1.1 --- a/benchmark/hcube.cc Wed Nov 14 17:53:08 2007 +0000
1.2 +++ b/benchmark/hcube.cc Sat Nov 17 20:58:11 2007 +0000
1.3 @@ -112,9 +112,9 @@
1.4 {
1.5 Graph::EdgeMap<int> flow(G);
1.6
1.7 - Preflow<Graph,int> MF(G,nodes[0],nodes[1<<dim-1],map,flow);
1.8 + Preflow<Graph> MF(G,map,nodes[0],nodes[1<<dim-1]);
1.9 for(int i=0;i<10;i++)
1.10 - MF.run(MF.NO_FLOW);
1.11 + MF.run();
1.12 }
1.13 PrintTime("PREFLOW",T);
1.14
2.1 --- a/demo/disjoint_paths_demo.cc Wed Nov 14 17:53:08 2007 +0000
2.2 +++ b/demo/disjoint_paths_demo.cc Sat Nov 17 20:58:11 2007 +0000
2.3 @@ -63,7 +63,7 @@
2.4
2.5 Graph::EdgeMap<int> flow(graph);
2.6
2.7 - Preflow<Graph, int, Capacity> preflow(graph, source, target, capacity, flow);
2.8 + Preflow<Graph, Capacity> preflow(graph, capacity, source, target);
2.9
2.10 preflow.run();
2.11
2.12 @@ -72,7 +72,7 @@
2.13 graphToEps(graph, "edge_disjoint_paths.eps").
2.14 title("edge disjoint paths").scaleToA4().
2.15 copyright("(C) 2003-2007 LEMON Project").drawArrows().
2.16 - edgeColors(composeMap(functorMap(color), flow)).
2.17 + edgeColors(composeMap(functorMap(color), preflow.flowMap())).
2.18 coords(coords).run();
2.19
2.20
2.21 @@ -87,9 +87,9 @@
2.22
2.23 SGraph::EdgeMap<int> sflow(sgraph);
2.24
2.25 - Preflow<SGraph, int, SCapacity> spreflow(sgraph, SGraph::outNode(source),
2.26 - SGraph::inNode(target),
2.27 - scapacity, sflow);
2.28 + Preflow<SGraph, SCapacity> spreflow(sgraph, scapacity,
2.29 + SGraph::outNode(source),
2.30 + SGraph::inNode(target));
2.31
2.32 spreflow.run();
2.33
2.34 @@ -99,7 +99,7 @@
2.35 graphToEps(sgraph, "node_disjoint_paths.eps").
2.36 title("node disjoint paths").scaleToA4().
2.37 copyright("(C) 2003-2007 LEMON Project").drawArrows().
2.38 - edgeColors(composeMap(functorMap(color), sflow)).
2.39 + edgeColors(composeMap(functorMap(color), spreflow.flowMap())).
2.40 coords(SGraph::combinedNodeMap(coords,
2.41 shiftMap(coords,
2.42 dim2::Point<double>(5, 0)))).
3.1 --- a/demo/sub_graph_adaptor_demo.cc Wed Nov 14 17:53:08 2007 +0000
3.2 +++ b/demo/sub_graph_adaptor_demo.cc Sat Nov 17 20:58:11 2007 +0000
3.3 @@ -109,10 +109,8 @@
3.4 SubGW gw(g, tight_edge_filter);
3.5
3.6 ConstMap<Edge, int> const_1_map(1);
3.7 - Graph::EdgeMap<int> flow(g, 0);
3.8 // Max flow between s and t in the graph of tight edges.
3.9 - Preflow<SubGW, int, ConstMap<Edge, int>, Graph::EdgeMap<int> >
3.10 - preflow(gw, s, t, const_1_map, flow);
3.11 + Preflow<SubGW, ConstMap<Edge, int> > preflow(gw, const_1_map, s, t);
3.12 preflow.run();
3.13
3.14 cout << "maximum number of edge-disjoint shortest paths: "
3.15 @@ -120,7 +118,7 @@
3.16 cout << "edges of the maximum number of edge-disjoint shortest s-t paths: "
3.17 << endl;
3.18 for(EdgeIt e(g); e!=INVALID; ++e)
3.19 - if (flow[e])
3.20 + if (preflow.flow(e) != 0)
3.21 cout << " " << g.id(e) << ", "
3.22 << g.id(g.source(e)) << "--"
3.23 << length[e] << "->" << g.id(g.target(e)) << endl;
4.1 --- a/doc/groups.dox Wed Nov 14 17:53:08 2007 +0000
4.2 +++ b/doc/groups.dox Sat Nov 17 20:58:11 2007 +0000
4.3 @@ -223,8 +223,27 @@
4.4 This group describes the algorithms for finding maximum flows and
4.5 feasible circulations.
4.6
4.7 -\image html flow.png
4.8 -\image latex flow.eps "Graph flow" width=\textwidth
4.9 +The maximum flow problem is to find a flow between a single-source and
4.10 +single-target that is maximum. Formally, there is \f$G=(V,A)\f$
4.11 +directed graph, an \f$c_a:A\rightarrow\mathbf{R}^+_0\f$ capacity
4.12 +function and given \f$s, t \in V\f$ source and target node. The
4.13 +maximum flow is the solution of the next optimization problem:
4.14 +
4.15 +\f[ 0 \le f_a \le c_a \f]
4.16 +\f[ \sum_{v\in\delta^{-}(u)}f_{vu}=\sum_{v\in\delta^{+}(u)}f_{uv} \quad u \in V \setminus \{s,t\}\f]
4.17 +\f[ \max \sum_{v\in\delta^{+}(s)}f_{uv} - \sum_{v\in\delta^{-}(s)}f_{vu}\f]
4.18 +
4.19 +The lemon contains several algorithms for solve maximum flow problems:
4.20 +- \ref lemon::EdmondsKarp "Edmonds-Karp"
4.21 +- \ref lemon::Preflow "Goldberg's Preflow algorithm"
4.22 +- \ref lemon::DinitzSleatorTarjan "Dinitz's blocking flow algorithm with dynamic tree"
4.23 +- \ref lemon::GoldbergTarjan "Preflow algorithm with dynamic trees"
4.24 +
4.25 +In most cases the \ref lemon::Preflow "preflow" algorithm provides the
4.26 +fastest method to compute the maximum flow. All impelementations
4.27 +provides functions for query the minimum cut, which is the dual linear
4.28 +programming probelm of the maximum flow.
4.29 +
4.30 */
4.31
4.32 /**
5.1 --- a/doc/images/flow.eps Wed Nov 14 17:53:08 2007 +0000
5.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
5.3 @@ -1,847 +0,0 @@
5.4 -%!PS-Adobe-3.0 EPSF-3.0
5.5 -%%BoundingBox: 13 59 827 531
5.6 -%%HiResBoundingBox: 13.1911 59.6331 826.078 530.254
5.7 -%%Creator: Karbon14 EPS Exportfilter 0.5
5.8 -%%CreationDate: (04/18/06 08:36:16)
5.9 -%%For: (Balazs Dezso) ()
5.10 -%%Title: ()
5.11 -
5.12 -/N {newpath} def
5.13 -/C {closepath} def
5.14 -/m {moveto} def
5.15 -/c {curveto} def
5.16 -/l {lineto} def
5.17 -/s {stroke} def
5.18 -/f {fill} def
5.19 -/w {setlinewidth} def
5.20 -/d {setdash} def
5.21 -/r {setrgbcolor} def
5.22 -/S {gsave} def
5.23 -/R {grestore} def
5.24 -
5.25 -N
5.26 -694.242 342.699 m
5.27 -748.648 255.656 753.472 247.937 799.086 174.957 c
5.28 -[] 0 d 0 0 0.5 r 3.58125 w s
5.29 -
5.30 -N
5.31 -807.312 161.797 m
5.32 -794.531 172.109 l
5.33 -803.64 177.805 l
5.34 -807.312 161.797 l
5.35 -C
5.36 -S 0 0 0.5 r f R
5.37 -
5.38 -N
5.39 -678.324 199.43 m
5.40 -736.031 179.062 744.519 176.066 787.746 160.812 c
5.41 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.42 -
5.43 -N
5.44 -802.382 155.644 m
5.45 -785.957 155.746 l
5.46 -789.531 165.879 l
5.47 -802.382 155.644 l
5.48 -C
5.49 -S 0.6 0 0.2 r f R
5.50 -
5.51 -N
5.52 -654.445 135.75 m
5.53 -723.707 142.676 732.703 143.578 786.316 148.937 c
5.54 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.55 -
5.56 -N
5.57 -801.761 150.484 m
5.58 -786.851 143.594 l
5.59 -785.781 154.285 l
5.60 -801.761 150.484 l
5.61 -C
5.62 -S 0.6 0 0.2 r f R
5.63 -
5.64 -N
5.65 -606.687 72.0741 m
5.66 -700.828 108.281 709.32 111.547 788.007 141.812 c
5.67 -[] 0 d 0.8 0 0.1 r 3.58125 w s
5.68 -
5.69 -N
5.70 -802.496 147.387 m
5.71 -789.937 136.797 l
5.72 -786.082 146.828 l
5.73 -802.496 147.387 l
5.74 -C
5.75 -S 0.8 0 0.1 r f R
5.76 -
5.77 -N
5.78 -694.242 342.699 m
5.79 -643.015 305.133 635.746 299.805 596.992 271.383 c
5.80 -[] 0 d 1 0 0 r 3.58125 w s
5.81 -
5.82 -N
5.83 -584.476 262.207 m
5.84 -593.816 275.715 l
5.85 -600.168 267.051 l
5.86 -584.476 262.207 l
5.87 -C
5.88 -S 1 0 0 r f R
5.89 -
5.90 -N
5.91 -582.808 358.621 m
5.92 -627.601 352.223 636.382 350.965 667.058 346.586 c
5.93 -[] 0 d 0.4 0 0.3 r 3.58125 w s
5.94 -
5.95 -N
5.96 -682.425 344.39 m
5.97 -666.3 341.266 l
5.98 -667.82 351.902 l
5.99 -682.425 344.39 l
5.100 -C
5.101 -S 0.4 0 0.3 r f R
5.102 -
5.103 -N
5.104 -511.171 517.816 m
5.105 -595.562 437.09 602.148 430.793 674.398 361.683 c
5.106 -[] 0 d 0.6 1 0.2 r 3.58125 w s
5.107 -
5.108 -N
5.109 -685.617 350.953 m
5.110 -670.687 357.801 l
5.111 -678.113 365.566 l
5.112 -685.617 350.953 l
5.113 -C
5.114 -S 0.6 1 0.2 r f R
5.115 -
5.116 -N
5.117 -495.25 143.711 m
5.118 -577.074 168.613 585.757 171.258 652.054 191.433 c
5.119 -[] 0 d 1 0 0 r 3.58125 w s
5.120 -
5.121 -N
5.122 -666.902 195.953 m
5.123 -653.617 186.293 l
5.124 -650.488 196.574 l
5.125 -666.902 195.953 l
5.126 -C
5.127 -S 1 0 0 r f R
5.128 -
5.129 -N
5.130 -574.847 255.148 m
5.131 -616.957 232.473 624.793 228.254 654.144 212.449 c
5.132 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.133 -
5.134 -N
5.135 -667.812 205.09 m
5.136 -651.597 207.719 l
5.137 -656.691 217.18 l
5.138 -667.812 205.09 l
5.139 -C
5.140 -S 0.6 0 0.2 r f R
5.141 -
5.142 -N
5.143 -495.25 143.711 m
5.144 -564.468 140.25 573.496 139.797 627.019 137.125 c
5.145 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.146 -
5.147 -N
5.148 -642.519 136.348 m
5.149 -626.75 131.758 l
5.150 -627.285 142.488 l
5.151 -642.519 136.348 l
5.152 -C
5.153 -S 0.6 0 0.2 r f R
5.154 -
5.155 -N
5.156 -391.773 72.0741 m
5.157 -489.199 72.0741 498.293 72.0741 579.226 72.0741 c
5.158 -[] 0 d 0.8 0 0.1 r 3.58125 w s
5.159 -
5.160 -N
5.161 -594.746 72.0741 m
5.162 -579.226 66.6991 l
5.163 -579.226 77.4451 l
5.164 -594.746 72.0741 l
5.165 -C
5.166 -S 0.8 0 0.1 r f R
5.167 -
5.168 -N
5.169 -495.25 143.711 m
5.170 -528.867 190.777 534.089 198.086 558.886 232.801 c
5.171 -[] 0 d 1 0 0 r 3.58125 w s
5.172 -
5.173 -N
5.174 -567.91 245.433 m
5.175 -563.257 229.68 l
5.176 -554.515 235.926 l
5.177 -567.91 245.433 l
5.178 -C
5.179 -S 1 0 0 r f R
5.180 -
5.181 -N
5.182 -431.574 271.062 m
5.183 -458.687 216.84 462.711 208.793 482.968 168.273 c
5.184 -[] 0 d 1 0 0 r 3.58125 w s
5.185 -
5.186 -N
5.187 -489.91 154.39 m
5.188 -478.164 165.871 l
5.189 -487.777 170.676 l
5.190 -489.91 154.39 l
5.191 -C
5.192 -S 1 0 0 r f R
5.193 -
5.194 -N
5.195 -391.773 72.0741 m
5.196 -434.636 101.75 441.992 106.84 472.671 128.082 c
5.197 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.198 -
5.199 -N
5.200 -485.433 136.914 m
5.201 -475.73 123.664 l
5.202 -469.613 132.496 l
5.203 -485.433 136.914 l
5.204 -C
5.205 -S 0.6 0 0.2 r f R
5.206 -
5.207 -N
5.208 -574.847 255.148 m
5.209 -577.964 295.68 578.64 304.457 580.703 331.238 c
5.210 -[] 0 d 1 0 0 r 3.58125 w s
5.211 -
5.212 -N
5.213 -581.89 346.715 m
5.214 -586.058 330.828 l
5.215 -575.343 331.652 l
5.216 -581.89 346.715 l
5.217 -C
5.218 -S 1 0 0 r f R
5.219 -
5.220 -N
5.221 -431.574 271.062 m
5.222 -492.73 264.269 501.675 263.273 547.554 258.18 c
5.223 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.224 -
5.225 -N
5.226 -562.98 256.465 m
5.227 -546.961 252.84 l
5.228 -548.148 263.516 l
5.229 -562.98 256.465 l
5.230 -C
5.231 -S 0.6 0 0.2 r f R
5.232 -
5.233 -N
5.234 -511.171 517.816 m
5.235 -542.777 447.578 546.492 439.32 571.539 383.664 c
5.236 -[] 0 d 0.6 1 0.2 r 3.58125 w s
5.237 -
5.238 -N
5.239 -577.906 369.508 m
5.240 -566.64 381.457 l
5.241 -576.437 385.867 l
5.242 -577.906 369.508 l
5.243 -C
5.244 -S 0.6 1 0.2 r f R
5.245 -
5.246 -N
5.247 -391.773 374.539 m
5.248 -477.171 367.422 486.218 366.668 555.441 360.902 c
5.249 -[] 0 d 0.8 1 0.1 r 3.58125 w s
5.250 -
5.251 -N
5.252 -570.91 359.613 m
5.253 -554.996 355.547 l
5.254 -555.886 366.254 l
5.255 -570.91 359.613 l
5.256 -C
5.257 -S 0.8 1 0.1 r f R
5.258 -
5.259 -N
5.260 -431.574 509.855 m
5.261 -459.371 512.637 467.753 513.473 483.843 515.082 c
5.262 -[] 0 d 0.2 0 0.4 r 3.58125 w s
5.263 -
5.264 -N
5.265 -499.289 516.625 m
5.266 -484.379 509.734 l
5.267 -483.312 520.43 l
5.268 -499.289 516.625 l
5.269 -C
5.270 -S 0.2 0 0.4 r f R
5.271 -
5.272 -N
5.273 -391.773 374.539 m
5.274 -444.953 438.351 450.761 445.324 493.589 496.719 c
5.275 -[] 0 d 1 0 0 r 3.58125 w s
5.276 -
5.277 -N
5.278 -503.527 508.645 m
5.279 -497.718 493.277 l
5.280 -489.461 500.16 l
5.281 -503.527 508.645 l
5.282 -C
5.283 -S 1 0 0 r f R
5.284 -
5.285 -N
5.286 -391.773 374.539 m
5.287 -408.687 432.043 411.226 440.676 423.824 483.508 c
5.288 -[] 0 d 1 0 0 r 3.58125 w s
5.289 -
5.290 -N
5.291 -428.203 498.402 m
5.292 -428.98 481.992 l
5.293 -418.668 485.027 l
5.294 -428.203 498.402 l
5.295 -C
5.296 -S 1 0 0 r f R
5.297 -
5.298 -N
5.299 -168.906 454.137 m
5.300 -290.617 479.953 299.535 481.848 404.711 504.156 c
5.301 -[] 0 d 0.2 0 0.4 r 3.58125 w s
5.302 -
5.303 -N
5.304 -419.894 507.379 m
5.305 -405.824 498.902 l
5.306 -403.593 509.414 l
5.307 -419.894 507.379 l
5.308 -C
5.309 -S 0.2 0 0.4 r f R
5.310 -
5.311 -N
5.312 -391.773 374.539 m
5.313 -407.699 333.133 410.879 324.863 421.714 296.695 c
5.314 -[] 0 d 0.6 1 0.2 r 3.58125 w s
5.315 -
5.316 -N
5.317 -427.285 282.207 m
5.318 -416.699 294.766 l
5.319 -426.73 298.621 l
5.320 -427.285 282.207 l
5.321 -C
5.322 -S 0.6 1 0.2 r f R
5.323 -
5.324 -N
5.325 -168.906 454.137 m
5.326 -270.98 417.683 279.554 414.621 365.918 383.777 c
5.327 -[] 0 d 0.4 0 0.3 r 3.58125 w s
5.328 -
5.329 -N
5.330 -380.535 378.555 m
5.331 -364.109 378.715 l
5.332 -367.722 388.836 l
5.333 -380.535 378.555 l
5.334 -C
5.335 -S 0.4 0 0.3 r f R
5.336 -
5.337 -N
5.338 -431.574 271.062 m
5.339 -392.125 238.195 385.261 232.473 357.156 209.051 c
5.340 -[] 0 d 1 1 0 r 3.58125 w s
5.341 -
5.342 -N
5.343 -345.23 199.113 m
5.344 -353.714 213.176 l
5.345 -360.597 204.922 l
5.346 -345.23 199.113 l
5.347 -C
5.348 -S 1 1 0 r f R
5.349 -
5.350 -N
5.351 -336.058 191.469 m
5.352 -359.39 141.473 363.183 133.348 380.164 96.9571 c
5.353 -[] 0 d 0.8 1 0.1 r 3.58125 w s
5.354 -
5.355 -N
5.356 -386.726 82.8941 m
5.357 -375.296 94.6871 l
5.358 -385.031 99.2301 l
5.359 -386.726 82.8941 l
5.360 -C
5.361 -S 0.8 1 0.1 r f R
5.362 -
5.363 -N
5.364 -208.703 143.711 m
5.365 -290.812 111.582 299.269 108.273 366.203 82.0821 c
5.366 -[] 0 d 0.6 1 0.2 r 3.58125 w s
5.367 -
5.368 -N
5.369 -380.66 76.4261 m
5.370 -364.246 77.0781 l
5.371 -368.164 87.0861 l
5.372 -380.66 76.4261 l
5.373 -C
5.374 -S 0.6 1 0.2 r f R
5.375 -
5.376 -N
5.377 -336.058 191.469 m
5.378 -297.203 262.117 292.836 270.055 261.738 326.598 c
5.379 -[] 0 d 0.8 0 0.1 r 3.58125 w s
5.380 -
5.381 -N
5.382 -254.257 340.199 m
5.383 -266.445 329.187 l
5.384 -257.031 324.008 l
5.385 -254.257 340.199 l
5.386 -C
5.387 -S 0.8 0 0.1 r f R
5.388 -
5.389 -N
5.390 -208.703 143.711 m
5.391 -262.414 163.851 270.82 167.004 310.347 181.828 c
5.392 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.393 -
5.394 -N
5.395 -324.879 187.277 m
5.396 -312.234 176.797 l
5.397 -308.461 186.859 l
5.398 -324.879 187.277 l
5.399 -C
5.400 -S 0.6 0 0.2 r f R
5.401 -
5.402 -N
5.403 -248.503 350.66 m
5.404 -215.242 393.902 209.777 401.004 185.648 432.371 c
5.405 -[] 0 d 0.4 0 0.3 r 3.58125 w s
5.406 -
5.407 -N
5.408 -176.183 444.672 m
5.409 -189.906 435.644 l
5.410 -181.39 429.094 l
5.411 -176.183 444.672 l
5.412 -C
5.413 -S 0.4 0 0.3 r f R
5.414 -
5.415 -N
5.416 -25.6286 271.062 m
5.417 -127.703 307.519 136.281 310.582 222.64 341.426 c
5.418 -[] 0 d 0.6 0 0.2 r 3.58125 w s
5.419 -
5.420 -N
5.421 -237.257 346.644 m
5.422 -224.449 336.363 l
5.423 -220.836 346.484 l
5.424 -237.257 346.644 l
5.425 -C
5.426 -S 0.6 0 0.2 r f R
5.427 -
5.428 -N
5.429 -25.6286 271.062 m
5.430 -91.1325 354.758 96.7419 361.93 151.98 432.512 c
5.431 -[] 0 d 0.2 0 0.4 r 3.58125 w s
5.432 -
5.433 -N
5.434 -161.546 444.734 m
5.435 -156.211 429.199 l
5.436 -147.75 435.824 l
5.437 -161.546 444.734 l
5.438 -C
5.439 -S 0.2 0 0.4 r f R
5.440 -
5.441 -N
5.442 -25.6286 271.062 m
5.443 -108.961 213.098 116.433 207.902 186.16 159.394 c
5.444 -[] 0 d 0.2 0 0.4 r 3.58125 w s
5.445 -
5.446 -N
5.447 -198.902 150.531 m
5.448 -183.093 154.984 l
5.449 -189.23 163.805 l
5.450 -198.902 150.531 l
5.451 -C
5.452 -S 0.2 0 0.4 r f R
5.453 -
5.454 -N
5.455 -37.57 271.062 m
5.456 -37.57 271.062 l
5.457 -37.57 277.656 32.2224 283.004 25.6286 283.004 c
5.458 -19.0348 283.004 13.6911 277.656 13.6911 271.062 c
5.459 -13.6911 264.469 19.0348 259.129 25.6286 259.129 c
5.460 -32.2224 259.129 37.57 264.469 37.57 271.062 c
5.461 -37.57 271.062 l
5.462 -C
5.463 -S 0 0 0 r f R
5.464 -
5.465 -N
5.466 -36.4841 271.062 m
5.467 -36.4841 271.062 l
5.468 -36.4841 277.055 31.6208 281.918 25.6286 281.918 c
5.469 -19.6364 281.918 14.777 277.055 14.777 271.062 c
5.470 -14.777 265.07 19.6364 260.215 25.6286 260.215 c
5.471 -31.6208 260.215 36.4841 265.07 36.4841 271.062 c
5.472 -36.4841 271.062 l
5.473 -C
5.474 -S 0 0 1 r f R
5.475 -
5.476 -N
5.477 -220.644 143.711 m
5.478 -220.644 143.711 l
5.479 -220.644 150.305 215.296 155.652 208.703 155.652 c
5.480 -202.109 155.652 196.765 150.305 196.765 143.711 c
5.481 -196.765 137.117 202.109 131.773 208.703 131.773 c
5.482 -215.296 131.773 220.644 137.117 220.644 143.711 c
5.483 -220.644 143.711 l
5.484 -C
5.485 -S 0 0 0 r f R
5.486 -
5.487 -N
5.488 -219.558 143.711 m
5.489 -219.558 143.711 l
5.490 -219.558 149.703 214.695 154.566 208.703 154.566 c
5.491 -202.711 154.566 197.851 149.703 197.851 143.711 c
5.492 -197.851 137.719 202.711 132.855 208.703 132.855 c
5.493 -214.695 132.855 219.558 137.719 219.558 143.711 c
5.494 -219.558 143.711 l
5.495 -C
5.496 -S 0.5 0.5 1 r f R
5.497 -
5.498 -N
5.499 -180.843 454.137 m
5.500 -180.843 454.137 l
5.501 -180.843 460.73 175.5 466.078 168.906 466.078 c
5.502 -162.312 466.078 156.964 460.73 156.964 454.137 c
5.503 -156.964 447.543 162.312 442.199 168.906 442.199 c
5.504 -175.5 442.199 180.843 447.543 180.843 454.137 c
5.505 -180.843 454.137 l
5.506 -C
5.507 -S 0 0 0 r f R
5.508 -
5.509 -N
5.510 -179.757 454.137 m
5.511 -179.757 454.137 l
5.512 -179.757 460.129 174.898 464.992 168.906 464.992 c
5.513 -162.914 464.992 158.05 460.129 158.05 454.137 c
5.514 -158.05 448.144 162.914 443.281 168.906 443.281 c
5.515 -174.898 443.281 179.757 448.144 179.757 454.137 c
5.516 -179.757 454.137 l
5.517 -C
5.518 -S 0.5 0.5 1 r f R
5.519 -
5.520 -N
5.521 -260.441 350.66 m
5.522 -260.441 350.66 l
5.523 -260.441 357.254 255.097 362.601 248.503 362.601 c
5.524 -241.91 362.601 236.562 357.254 236.562 350.66 c
5.525 -236.562 344.066 241.91 338.723 248.503 338.723 c
5.526 -255.097 338.723 260.441 344.066 260.441 350.66 c
5.527 -260.441 350.66 l
5.528 -C
5.529 -S 0 0 0 r f R
5.530 -
5.531 -N
5.532 -259.359 350.66 m
5.533 -259.359 350.66 l
5.534 -259.359 356.652 254.496 361.516 248.503 361.516 c
5.535 -242.511 361.516 237.648 356.652 237.648 350.66 c
5.536 -237.648 344.668 242.511 339.805 248.503 339.805 c
5.537 -254.496 339.805 259.359 344.668 259.359 350.66 c
5.538 -259.359 350.66 l
5.539 -C
5.540 -S 0.5 0.5 1 r f R
5.541 -
5.542 -N
5.543 -348 191.469 m
5.544 -348 191.469 l
5.545 -348 198.062 342.652 203.41 336.058 203.41 c
5.546 -329.464 203.41 324.121 198.062 324.121 191.469 c
5.547 -324.121 184.875 329.464 179.531 336.058 179.531 c
5.548 -342.652 179.531 348 184.875 348 191.469 c
5.549 -348 191.469 l
5.550 -C
5.551 -S 0 0 0 r f R
5.552 -
5.553 -N
5.554 -346.914 191.469 m
5.555 -346.914 191.469 l
5.556 -346.914 197.461 342.05 202.324 336.058 202.324 c
5.557 -330.066 202.324 325.207 197.461 325.207 191.469 c
5.558 -325.207 185.476 330.066 180.617 336.058 180.617 c
5.559 -342.05 180.617 346.914 185.476 346.914 191.469 c
5.560 -346.914 191.469 l
5.561 -C
5.562 -S 0.5 0.5 1 r f R
5.563 -
5.564 -N
5.565 -403.714 72.0741 m
5.566 -403.714 72.0741 l
5.567 -403.714 78.6681 398.367 84.0121 391.773 84.0121 c
5.568 -385.179 84.0121 379.839 78.6681 379.839 72.0741 c
5.569 -379.839 65.4801 385.179 60.1331 391.773 60.1331 c
5.570 -398.367 60.1331 403.714 65.4801 403.714 72.0741 c
5.571 -403.714 72.0741 l
5.572 -C
5.573 -S 0 0 0 r f R
5.574 -
5.575 -N
5.576 -402.629 72.0741 m
5.577 -402.629 72.0741 l
5.578 -402.629 78.0661 397.765 82.9301 391.773 82.9301 c
5.579 -385.781 82.9301 380.925 78.0661 380.925 72.0741 c
5.580 -380.925 66.0821 385.781 61.2191 391.773 61.2191 c
5.581 -397.765 61.2191 402.629 66.0821 402.629 72.0741 c
5.582 -402.629 72.0741 l
5.583 -C
5.584 -S 1 0.5 1 r f R
5.585 -
5.586 -N
5.587 -443.511 271.062 m
5.588 -443.511 271.062 l
5.589 -443.511 277.656 438.168 283.004 431.574 283.004 c
5.590 -424.98 283.004 419.632 277.656 419.632 271.062 c
5.591 -419.632 264.469 424.98 259.129 431.574 259.129 c
5.592 -438.168 259.129 443.511 264.469 443.511 271.062 c
5.593 -443.511 271.062 l
5.594 -C
5.595 -S 0 0 0 r f R
5.596 -
5.597 -N
5.598 -442.425 271.062 m
5.599 -442.425 271.062 l
5.600 -442.425 277.055 437.566 281.918 431.574 281.918 c
5.601 -425.582 281.918 420.718 277.055 420.718 271.062 c
5.602 -420.718 265.07 425.582 260.215 431.574 260.215 c
5.603 -437.566 260.215 442.425 265.07 442.425 271.062 c
5.604 -442.425 271.062 l
5.605 -C
5.606 -S 1 0.5 1 r f R
5.607 -
5.608 -N
5.609 -403.714 374.539 m
5.610 -403.714 374.539 l
5.611 -403.714 381.133 398.367 386.48 391.773 386.48 c
5.612 -385.179 386.48 379.839 381.133 379.839 374.539 c
5.613 -379.839 367.945 385.179 362.601 391.773 362.601 c
5.614 -398.367 362.601 403.714 367.945 403.714 374.539 c
5.615 -403.714 374.539 l
5.616 -C
5.617 -S 0 0 0 r f R
5.618 -
5.619 -N
5.620 -402.629 374.539 m
5.621 -402.629 374.539 l
5.622 -402.629 380.531 397.765 385.394 391.773 385.394 c
5.623 -385.781 385.394 380.925 380.531 380.925 374.539 c
5.624 -380.925 368.547 385.781 363.687 391.773 363.687 c
5.625 -397.765 363.687 402.629 368.547 402.629 374.539 c
5.626 -402.629 374.539 l
5.627 -C
5.628 -S 0.5 0.5 1 r f R
5.629 -
5.630 -N
5.631 -443.511 509.855 m
5.632 -443.511 509.855 l
5.633 -443.511 516.449 438.168 521.797 431.574 521.797 c
5.634 -424.98 521.797 419.632 516.449 419.632 509.855 c
5.635 -419.632 503.262 424.98 497.918 431.574 497.918 c
5.636 -438.168 497.918 443.511 503.262 443.511 509.855 c
5.637 -443.511 509.855 l
5.638 -C
5.639 -S 0 0 0 r f R
5.640 -
5.641 -N
5.642 -442.425 509.855 m
5.643 -442.425 509.855 l
5.644 -442.425 515.848 437.566 520.711 431.574 520.711 c
5.645 -425.582 520.711 420.718 515.848 420.718 509.855 c
5.646 -420.718 503.863 425.582 499 431.574 499 c
5.647 -437.566 499 442.425 503.863 442.425 509.855 c
5.648 -442.425 509.855 l
5.649 -C
5.650 -S 0.5 0.5 1 r f R
5.651 -
5.652 -N
5.653 -523.109 517.816 m
5.654 -523.109 517.816 l
5.655 -523.109 524.41 517.765 529.754 511.171 529.754 c
5.656 -504.578 529.754 499.23 524.41 499.23 517.816 c
5.657 -499.23 511.223 504.578 505.875 511.171 505.875 c
5.658 -517.765 505.875 523.109 511.223 523.109 517.816 c
5.659 -523.109 517.816 l
5.660 -C
5.661 -S 0 0 0 r f R
5.662 -
5.663 -N
5.664 -522.023 517.816 m
5.665 -522.023 517.816 l
5.666 -522.023 523.809 517.164 528.668 511.171 528.668 c
5.667 -505.179 528.668 500.316 523.809 500.316 517.816 c
5.668 -500.316 511.824 505.179 506.961 511.171 506.961 c
5.669 -517.164 506.961 522.023 511.824 522.023 517.816 c
5.670 -522.023 517.816 l
5.671 -C
5.672 -S 0.5 0.5 1 r f R
5.673 -
5.674 -N
5.675 -594.746 358.621 m
5.676 -594.746 358.621 l
5.677 -594.746 365.215 589.402 370.558 582.808 370.558 c
5.678 -576.214 370.558 570.867 365.215 570.867 358.621 c
5.679 -570.867 352.027 576.214 346.68 582.808 346.68 c
5.680 -589.402 346.68 594.746 352.027 594.746 358.621 c
5.681 -594.746 358.621 l
5.682 -C
5.683 -S 0 0 0 r f R
5.684 -
5.685 -N
5.686 -593.664 358.621 m
5.687 -593.664 358.621 l
5.688 -593.664 364.613 588.8 369.473 582.808 369.473 c
5.689 -576.816 369.473 571.953 364.613 571.953 358.621 c
5.690 -571.953 352.629 576.816 347.766 582.808 347.766 c
5.691 -588.8 347.766 593.664 352.629 593.664 358.621 c
5.692 -593.664 358.621 l
5.693 -C
5.694 -S 1 0.5 1 r f R
5.695 -
5.696 -N
5.697 -586.789 255.148 m
5.698 -586.789 255.148 l
5.699 -586.789 261.742 581.441 267.082 574.847 267.082 c
5.700 -568.253 267.082 562.91 261.742 562.91 255.148 c
5.701 -562.91 248.555 568.253 243.207 574.847 243.207 c
5.702 -581.441 243.207 586.789 248.555 586.789 255.148 c
5.703 -586.789 255.148 l
5.704 -C
5.705 -S 0 0 0 r f R
5.706 -
5.707 -N
5.708 -585.703 255.148 m
5.709 -585.703 255.148 l
5.710 -585.703 261.14 580.839 265.996 574.847 265.996 c
5.711 -568.855 265.996 563.992 261.14 563.992 255.148 c
5.712 -563.992 249.156 568.855 244.293 574.847 244.293 c
5.713 -580.839 244.293 585.703 249.156 585.703 255.148 c
5.714 -585.703 255.148 l
5.715 -C
5.716 -S 1 0.5 1 r f R
5.717 -
5.718 -N
5.719 -507.191 143.711 m
5.720 -507.191 143.711 l
5.721 -507.191 150.305 501.843 155.652 495.25 155.652 c
5.722 -488.656 155.652 483.312 150.305 483.312 143.711 c
5.723 -483.312 137.117 488.656 131.773 495.25 131.773 c
5.724 -501.843 131.773 507.191 137.117 507.191 143.711 c
5.725 -507.191 143.711 l
5.726 -C
5.727 -S 0 0 0 r f R
5.728 -
5.729 -N
5.730 -506.105 143.711 m
5.731 -506.105 143.711 l
5.732 -506.105 149.703 501.242 154.566 495.25 154.566 c
5.733 -489.257 154.566 484.398 149.703 484.398 143.711 c
5.734 -484.398 137.719 489.257 132.855 495.25 132.855 c
5.735 -501.242 132.855 506.105 137.719 506.105 143.711 c
5.736 -506.105 143.711 l
5.737 -C
5.738 -S 1 0.5 1 r f R
5.739 -
5.740 -N
5.741 -618.629 72.0741 m
5.742 -618.629 72.0741 l
5.743 -618.629 78.6681 613.281 84.0121 606.687 84.0121 c
5.744 -600.093 84.0121 594.746 78.6681 594.746 72.0741 c
5.745 -594.746 65.4801 600.093 60.1331 606.687 60.1331 c
5.746 -613.281 60.1331 618.629 65.4801 618.629 72.0741 c
5.747 -618.629 72.0741 l
5.748 -C
5.749 -S 0 0 0 r f R
5.750 -
5.751 -N
5.752 -617.543 72.0741 m
5.753 -617.543 72.0741 l
5.754 -617.543 78.0661 612.679 82.9301 606.687 82.9301 c
5.755 -600.695 82.9301 595.832 78.0661 595.832 72.0741 c
5.756 -595.832 66.0821 600.695 61.2191 606.687 61.2191 c
5.757 -612.679 61.2191 617.543 66.0821 617.543 72.0741 c
5.758 -617.543 72.0741 l
5.759 -C
5.760 -S 1 0.5 1 r f R
5.761 -
5.762 -N
5.763 -666.386 135.75 m
5.764 -666.386 135.75 l
5.765 -666.386 142.344 661.039 147.691 654.445 147.691 c
5.766 -647.851 147.691 642.507 142.344 642.507 135.75 c
5.767 -642.507 129.156 647.851 123.812 654.445 123.812 c
5.768 -661.039 123.812 666.386 129.156 666.386 135.75 c
5.769 -666.386 135.75 l
5.770 -C
5.771 -S 0 0 0 r f R
5.772 -
5.773 -N
5.774 -665.3 135.75 m
5.775 -665.3 135.75 l
5.776 -665.3 141.742 660.437 146.605 654.445 146.605 c
5.777 -648.453 146.605 643.589 141.742 643.589 135.75 c
5.778 -643.589 129.758 648.453 124.898 654.445 124.898 c
5.779 -660.437 124.898 665.3 129.758 665.3 135.75 c
5.780 -665.3 135.75 l
5.781 -C
5.782 -S 1 0.5 1 r f R
5.783 -
5.784 -N
5.785 -690.265 199.43 m
5.786 -690.265 199.43 l
5.787 -690.265 206.023 684.918 211.371 678.324 211.371 c
5.788 -671.73 211.371 666.386 206.023 666.386 199.43 c
5.789 -666.386 192.836 671.73 187.488 678.324 187.488 c
5.790 -684.918 187.488 690.265 192.836 690.265 199.43 c
5.791 -690.265 199.43 l
5.792 -C
5.793 -S 0 0 0 r f R
5.794 -
5.795 -N
5.796 -689.179 199.43 m
5.797 -689.179 199.43 l
5.798 -689.179 205.422 684.316 210.285 678.324 210.285 c
5.799 -672.332 210.285 667.472 205.422 667.472 199.43 c
5.800 -667.472 193.437 672.332 188.574 678.324 188.574 c
5.801 -684.316 188.574 689.179 193.437 689.179 199.43 c
5.802 -689.179 199.43 l
5.803 -C
5.804 -S 1 0.5 1 r f R
5.805 -
5.806 -N
5.807 -706.183 342.699 m
5.808 -706.183 342.699 l
5.809 -706.183 349.293 700.836 354.64 694.242 354.64 c
5.810 -687.648 354.64 682.304 349.293 682.304 342.699 c
5.811 -682.304 336.105 687.648 330.762 694.242 330.762 c
5.812 -700.836 330.762 706.183 336.105 706.183 342.699 c
5.813 -706.183 342.699 l
5.814 -C
5.815 -S 0 0 0 r f R
5.816 -
5.817 -N
5.818 -705.097 342.699 m
5.819 -705.097 342.699 l
5.820 -705.097 348.691 700.234 353.555 694.242 353.555 c
5.821 -688.25 353.555 683.39 348.691 683.39 342.699 c
5.822 -683.39 336.707 688.25 331.848 694.242 331.848 c
5.823 -700.234 331.848 705.097 336.707 705.097 342.699 c
5.824 -705.097 342.699 l
5.825 -C
5.826 -S 1 0.5 1 r f R
5.827 -
5.828 -N
5.829 -825.578 151.672 m
5.830 -825.578 151.672 l
5.831 -825.578 158.266 820.234 163.609 813.64 163.609 c
5.832 -807.046 163.609 801.699 158.266 801.699 151.672 c
5.833 -801.699 145.078 807.046 139.73 813.64 139.73 c
5.834 -820.234 139.73 825.578 145.078 825.578 151.672 c
5.835 -825.578 151.672 l
5.836 -C
5.837 -S 0 0 0 r f R
5.838 -
5.839 -N
5.840 -824.496 151.672 m
5.841 -824.496 151.672 l
5.842 -824.496 157.664 819.632 162.523 813.64 162.523 c
5.843 -807.648 162.523 802.785 157.664 802.785 151.672 c
5.844 -802.785 145.68 807.648 140.816 813.64 140.816 c
5.845 -819.632 140.816 824.496 145.68 824.496 151.672 c
5.846 -824.496 151.672 l
5.847 -C
5.848 -S 1 0 1 r f R
5.849 -
5.850 -%%EOF
6.1 Binary file doc/images/flow.png has changed
7.1 --- a/lemon/Makefile.am Wed Nov 14 17:53:08 2007 +0000
7.2 +++ b/lemon/Makefile.am Sat Nov 17 20:58:11 2007 +0000
7.3 @@ -52,9 +52,11 @@
7.4 lemon/dag_shortest_path.h \
7.5 lemon/dfs.h \
7.6 lemon/dijkstra.h \
7.7 + lemon/dinitz_sleator_tarjan.h \
7.8 lemon/dist_log.h \
7.9 lemon/dim2.h \
7.10 lemon/dimacs.h \
7.11 + lemon/dynamic_tree.h \
7.12 lemon/edge_set.h \
7.13 lemon/edmonds_karp.h \
7.14 lemon/elevator.h \
7.15 @@ -71,6 +73,7 @@
7.16 lemon/graph_utils.h \
7.17 lemon/graph_writer.h \
7.18 lemon/grid_ugraph.h \
7.19 + lemon/goldberg_tarjan.h \
7.20 lemon/hao_orlin.h \
7.21 lemon/hypercube_graph.h \
7.22 lemon/iterable_maps.h \
8.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
8.2 +++ b/lemon/dinitz_sleator_tarjan.h Sat Nov 17 20:58:11 2007 +0000
8.3 @@ -0,0 +1,762 @@
8.4 +/* -*- C++ -*-
8.5 + *
8.6 + * This file is a part of LEMON, a generic C++ optimization library
8.7 + *
8.8 + * Copyright (C) 2003-2007
8.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
8.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
8.11 + *
8.12 + * Permission to use, modify and distribute this software is granted
8.13 + * provided that this copyright notice appears in all copies. For
8.14 + * precise terms see the accompanying LICENSE file.
8.15 + *
8.16 + * This software is provided "AS IS" with no warranty of any kind,
8.17 + * express or implied, and with no claim as to its suitability for any
8.18 + * purpose.
8.19 + *
8.20 + */
8.21 +
8.22 +#ifndef LEMON_DINITZ_SLEATOR_TARJAN_H
8.23 +#define LEMON_DINITZ_SLEATOR_TARJAN_H
8.24 +
8.25 +/// \file
8.26 +/// \ingroup max_flow
8.27 +/// \brief Implementation the dynamic tree data structure of Sleator
8.28 +/// and Tarjan.
8.29 +
8.30 +#include <lemon/time_measure.h>
8.31 +#include <queue>
8.32 +#include <lemon/graph_utils.h>
8.33 +#include <lemon/tolerance.h>
8.34 +#include <lemon/graph_adaptor.h>
8.35 +#include <lemon/bfs.h>
8.36 +#include <lemon/edge_set.h>
8.37 +#include <lemon/dynamic_tree.h>
8.38 +
8.39 +#include <vector>
8.40 +#include <limits>
8.41 +#include <fstream>
8.42 +
8.43 +
8.44 +namespace lemon {
8.45 +
8.46 + /// \brief Default traits class of DinitzSleatorTarjan class.
8.47 + ///
8.48 + /// Default traits class of DinitzSleatorTarjan class.
8.49 + /// \param _Graph Graph type.
8.50 + /// \param _CapacityMap Type of capacity map.
8.51 + template <typename _Graph, typename _CapacityMap>
8.52 + struct DinitzSleatorTarjanDefaultTraits {
8.53 +
8.54 + /// \brief The graph type the algorithm runs on.
8.55 + typedef _Graph Graph;
8.56 +
8.57 + /// \brief The type of the map that stores the edge capacities.
8.58 + ///
8.59 + /// The type of the map that stores the edge capacities.
8.60 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
8.61 + typedef _CapacityMap CapacityMap;
8.62 +
8.63 + /// \brief The type of the length of the edges.
8.64 + typedef typename CapacityMap::Value Value;
8.65 +
8.66 + /// \brief The map type that stores the flow values.
8.67 + ///
8.68 + /// The map type that stores the flow values.
8.69 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
8.70 + typedef typename Graph::template EdgeMap<Value> FlowMap;
8.71 +
8.72 + /// \brief Instantiates a FlowMap.
8.73 + ///
8.74 + /// This function instantiates a \ref FlowMap.
8.75 + /// \param graph The graph, to which we would like to define the flow map.
8.76 + static FlowMap* createFlowMap(const Graph& graph) {
8.77 + return new FlowMap(graph);
8.78 + }
8.79 +
8.80 + /// \brief The tolerance used by the algorithm
8.81 + ///
8.82 + /// The tolerance used by the algorithm to handle inexact computation.
8.83 + typedef Tolerance<Value> Tolerance;
8.84 +
8.85 + };
8.86 +
8.87 + /// \ingroup max_flow
8.88 + ///
8.89 + /// \brief Dinitz-Sleator-Tarjan algorithms class.
8.90 + ///
8.91 + /// This class provides an implementation of the \e
8.92 + /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
8.93 + /// value in a directed graph. The DinitzSleatorTarjan algorithm is
8.94 + /// the fastest known max flow algorithms wich using blocking
8.95 + /// flow. It is an improvement of the Dinitz algorithm by using the
8.96 + /// \ref DynamicTree "dynamic tree" data structure of Sleator and
8.97 + /// Tarjan.
8.98 + ///
8.99 + /// This blocking flow algorithms builds a layered graph according
8.100 + /// to \ref Bfs "breadth-first search" distance from the target node
8.101 + /// in the reversed residual graph. The layered graph contains each
8.102 + /// residual edge which steps one level down. After that the
8.103 + /// algorithm constructs a blocking flow on the layered graph and
8.104 + /// augments the overall flow with it. The number of the levels in
8.105 + /// the layered graph is strictly increasing in each augmenting
8.106 + /// phase therefore the number of the augmentings is at most
8.107 + /// \f$n-1\f$. The length of each phase is at most
8.108 + /// \f$O(m\log(n))\f$, that the overall time complexity is
8.109 + /// \f$O(nm\log(n))\f$.
8.110 + ///
8.111 + /// \param _Graph The directed graph type the algorithm runs on.
8.112 + /// \param _CapacityMap The capacity map type.
8.113 + /// \param _Traits Traits class to set various data types used by
8.114 + /// the algorithm. The default traits class is \ref
8.115 + /// DinitzSleatorTarjanDefaultTraits. See \ref
8.116 + /// DinitzSleatorTarjanDefaultTraits for the documentation of a
8.117 + /// Dinitz-Sleator-Tarjan traits class.
8.118 + ///
8.119 + /// \author Tamas Hamori and Balazs Dezso
8.120 +#ifdef DOXYGEN
8.121 + template <typename _Graph, typename _CapacityMap, typename _Traits>
8.122 +#else
8.123 + template <typename _Graph,
8.124 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
8.125 + typename _Traits =
8.126 + DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
8.127 +#endif
8.128 + class DinitzSleatorTarjan {
8.129 + public:
8.130 +
8.131 + typedef _Traits Traits;
8.132 + typedef typename Traits::Graph Graph;
8.133 + typedef typename Traits::CapacityMap CapacityMap;
8.134 + typedef typename Traits::Value Value;
8.135 +
8.136 + typedef typename Traits::FlowMap FlowMap;
8.137 + typedef typename Traits::Tolerance Tolerance;
8.138 +
8.139 +
8.140 + private:
8.141 +
8.142 + GRAPH_TYPEDEFS(typename Graph);
8.143 +
8.144 +
8.145 + typedef typename Graph::template NodeMap<int> LevelMap;
8.146 + typedef typename Graph::template NodeMap<int> IntNodeMap;
8.147 + typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
8.148 + typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
8.149 +
8.150 + private:
8.151 +
8.152 + const Graph& _graph;
8.153 + const CapacityMap* _capacity;
8.154 +
8.155 + Node _source, _target;
8.156 +
8.157 + FlowMap* _flow;
8.158 + bool _local_flow;
8.159 +
8.160 + IntNodeMap* _level;
8.161 + EdgeNodeMap* _dt_edges;
8.162 +
8.163 + IntNodeMap* _dt_index;
8.164 + DynTree* _dt;
8.165 +
8.166 + Tolerance _tolerance;
8.167 +
8.168 + Value _flow_value;
8.169 + Value _max_value;
8.170 +
8.171 +
8.172 + public:
8.173 +
8.174 + typedef DinitzSleatorTarjan Create;
8.175 +
8.176 + ///\name Named template parameters
8.177 +
8.178 + ///@{
8.179 +
8.180 + template <typename _FlowMap>
8.181 + struct DefFlowMapTraits : public Traits {
8.182 + typedef _FlowMap FlowMap;
8.183 + static FlowMap *createFlowMap(const Graph&) {
8.184 + throw UninitializedParameter();
8.185 + }
8.186 + };
8.187 +
8.188 + /// \brief \ref named-templ-param "Named parameter" for setting
8.189 + /// FlowMap type
8.190 + ///
8.191 + /// \ref named-templ-param "Named parameter" for setting FlowMap
8.192 + /// type
8.193 + template <typename _FlowMap>
8.194 + struct DefFlowMap
8.195 + : public DinitzSleatorTarjan<Graph, CapacityMap,
8.196 + DefFlowMapTraits<_FlowMap> > {
8.197 + typedef DinitzSleatorTarjan<Graph, CapacityMap,
8.198 + DefFlowMapTraits<_FlowMap> > Create;
8.199 + };
8.200 +
8.201 + template <typename _Elevator>
8.202 + struct DefElevatorTraits : public Traits {
8.203 + typedef _Elevator Elevator;
8.204 + static Elevator *createElevator(const Graph&, int) {
8.205 + throw UninitializedParameter();
8.206 + }
8.207 + };
8.208 +
8.209 + /// @}
8.210 +
8.211 + /// \brief \ref Exception for the case when the source equals the target.
8.212 + ///
8.213 + /// \ref Exception for the case when the source equals the target.
8.214 + ///
8.215 + class InvalidArgument : public lemon::LogicError {
8.216 + public:
8.217 + virtual const char* what() const throw() {
8.218 + return "lemon::DinitzSleatorTarjan::InvalidArgument";
8.219 + }
8.220 + };
8.221 +
8.222 + /// \brief The constructor of the class.
8.223 + ///
8.224 + /// The constructor of the class.
8.225 + /// \param graph The directed graph the algorithm runs on.
8.226 + /// \param capacity The capacity of the edges.
8.227 + /// \param source The source node.
8.228 + /// \param target The target node.
8.229 + DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
8.230 + Node source, Node target)
8.231 + : _graph(graph), _capacity(&capacity),
8.232 + _source(source), _target(target),
8.233 + _flow(0), _local_flow(false),
8.234 + _level(0), _dt_edges(0),
8.235 + _dt_index(0), _dt(0),
8.236 + _tolerance(), _flow_value(), _max_value()
8.237 + {
8.238 + if (_source == _target) {
8.239 + throw InvalidArgument();
8.240 + }
8.241 + }
8.242 +
8.243 + /// \brief Destrcutor.
8.244 + ///
8.245 + /// Destructor.
8.246 + ~DinitzSleatorTarjan() {
8.247 + destroyStructures();
8.248 + }
8.249 +
8.250 + /// \brief Sets the capacity map.
8.251 + ///
8.252 + /// Sets the capacity map.
8.253 + /// \return \c (*this)
8.254 + DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
8.255 + _capacity = ↦
8.256 + return *this;
8.257 + }
8.258 +
8.259 + /// \brief Sets the flow map.
8.260 + ///
8.261 + /// Sets the flow map.
8.262 + /// \return \c (*this)
8.263 + DinitzSleatorTarjan& flowMap(FlowMap& map) {
8.264 + if (_local_flow) {
8.265 + delete _flow;
8.266 + _local_flow = false;
8.267 + }
8.268 + _flow = ↦
8.269 + return *this;
8.270 + }
8.271 +
8.272 + /// \brief Returns the flow map.
8.273 + ///
8.274 + /// \return The flow map.
8.275 + const FlowMap& flowMap() {
8.276 + return *_flow;
8.277 + }
8.278 +
8.279 + /// \brief Sets the source node.
8.280 + ///
8.281 + /// Sets the source node.
8.282 + /// \return \c (*this)
8.283 + DinitzSleatorTarjan& source(const Node& node) {
8.284 + _source = node;
8.285 + return *this;
8.286 + }
8.287 +
8.288 + /// \brief Sets the target node.
8.289 + ///
8.290 + /// Sets the target node.
8.291 + /// \return \c (*this)
8.292 + DinitzSleatorTarjan& target(const Node& node) {
8.293 + _target = node;
8.294 + return *this;
8.295 + }
8.296 +
8.297 + /// \brief Sets the tolerance used by algorithm.
8.298 + ///
8.299 + /// Sets the tolerance used by algorithm.
8.300 + DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
8.301 + _tolerance = tolerance;
8.302 + if (_dt) {
8.303 + _dt.tolerance(_tolerance);
8.304 + }
8.305 + return *this;
8.306 + }
8.307 +
8.308 + /// \brief Returns the tolerance used by algorithm.
8.309 + ///
8.310 + /// Returns the tolerance used by algorithm.
8.311 + const Tolerance& tolerance() const {
8.312 + return tolerance;
8.313 + }
8.314 +
8.315 + private:
8.316 +
8.317 + void createStructures() {
8.318 + if (!_flow) {
8.319 + _flow = Traits::createFlowMap(_graph);
8.320 + _local_flow = true;
8.321 + }
8.322 + if (!_level) {
8.323 + _level = new LevelMap(_graph);
8.324 + }
8.325 + if (!_dt_index && !_dt) {
8.326 + _dt_index = new IntNodeMap(_graph);
8.327 + _dt = new DynTree(*_dt_index, _tolerance);
8.328 + }
8.329 + if (!_dt_edges) {
8.330 + _dt_edges = new EdgeNodeMap(_graph);
8.331 + }
8.332 + _max_value = _dt->maxValue();
8.333 + }
8.334 +
8.335 + void destroyStructures() {
8.336 + if (_local_flow) {
8.337 + delete _flow;
8.338 + }
8.339 + if (_level) {
8.340 + delete _level;
8.341 + }
8.342 + if (_dt) {
8.343 + delete _dt;
8.344 + }
8.345 + if (_dt_index) {
8.346 + delete _dt_index;
8.347 + }
8.348 + if (_dt_edges) {
8.349 + delete _dt_edges;
8.350 + }
8.351 + }
8.352 +
8.353 + bool createLayeredGraph() {
8.354 +
8.355 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.356 + _level->set(n, -2);
8.357 + }
8.358 +
8.359 + int level = 0;
8.360 +
8.361 + std::vector<Node> queue;
8.362 + queue.push_back(_target);
8.363 + _level->set(_target, level);
8.364 +
8.365 + while ((*_level)[_source] == -2 && !queue.empty()) {
8.366 + std::vector<Node> nqueue;
8.367 + ++level;
8.368 +
8.369 + for (int i = 0; i < int(queue.size()); ++i) {
8.370 + Node n = queue[i];
8.371 +
8.372 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
8.373 + Node v = _graph.target(e);
8.374 + if ((*_level)[v] != -2) continue;
8.375 + Value rem = (*_flow)[e];
8.376 + if (!_tolerance.positive(rem)) continue;
8.377 + _level->set(v, level);
8.378 + nqueue.push_back(v);
8.379 + }
8.380 +
8.381 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
8.382 + Node v = _graph.source(e);
8.383 + if ((*_level)[v] != -2) continue;
8.384 + Value rem = (*_capacity)[e] - (*_flow)[e];
8.385 + if (!_tolerance.positive(rem)) continue;
8.386 + _level->set(v, level);
8.387 + nqueue.push_back(v);
8.388 + }
8.389 +
8.390 + }
8.391 + queue.swap(nqueue);
8.392 + }
8.393 +
8.394 + return (*_level)[_source] != -2;
8.395 + }
8.396 +
8.397 + void initEdges() {
8.398 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.399 + _graph.firstOut((*_dt_edges)[n], n);
8.400 + }
8.401 + }
8.402 +
8.403 +
8.404 + void augmentPath() {
8.405 + Value rem;
8.406 + Node n = _dt->findCost(_source, rem);
8.407 + _flow_value += rem;
8.408 + _dt->addCost(_source, - rem);
8.409 +
8.410 + _dt->cut(n);
8.411 + _dt->addCost(n, _max_value);
8.412 +
8.413 + Edge e = (*_dt_edges)[n];
8.414 + if (_graph.source(e) == n) {
8.415 + _flow->set(e, (*_capacity)[e]);
8.416 +
8.417 + _graph.nextOut(e);
8.418 + if (e == INVALID) {
8.419 + _graph.firstIn(e, n);
8.420 + }
8.421 + } else {
8.422 + _flow->set(e, 0);
8.423 + _graph.nextIn(e);
8.424 + }
8.425 + _dt_edges->set(n, e);
8.426 +
8.427 + }
8.428 +
8.429 + bool advance(Node n) {
8.430 + Edge e = (*_dt_edges)[n];
8.431 + if (e == INVALID) return false;
8.432 +
8.433 + Node u;
8.434 + Value rem;
8.435 + if (_graph.source(e) == n) {
8.436 + u = _graph.target(e);
8.437 + while ((*_level)[n] != (*_level)[u] + 1 ||
8.438 + !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
8.439 + _graph.nextOut(e);
8.440 + if (e == INVALID) break;
8.441 + u = _graph.target(e);
8.442 + }
8.443 + if (e != INVALID) {
8.444 + rem = (*_capacity)[e] - (*_flow)[e];
8.445 + } else {
8.446 + _graph.firstIn(e, n);
8.447 + if (e == INVALID) {
8.448 + _dt_edges->set(n, INVALID);
8.449 + return false;
8.450 + }
8.451 + u = _graph.source(e);
8.452 + while ((*_level)[n] != (*_level)[u] + 1 ||
8.453 + !_tolerance.positive((*_flow)[e])) {
8.454 + _graph.nextIn(e);
8.455 + if (e == INVALID) {
8.456 + _dt_edges->set(n, INVALID);
8.457 + return false;
8.458 + }
8.459 + u = _graph.source(e);
8.460 + }
8.461 + rem = (*_flow)[e];
8.462 + }
8.463 + } else {
8.464 + u = _graph.source(e);
8.465 + while ((*_level)[n] != (*_level)[u] + 1 ||
8.466 + !_tolerance.positive((*_flow)[e])) {
8.467 + _graph.nextIn(e);
8.468 + if (e == INVALID) {
8.469 + _dt_edges->set(n, INVALID);
8.470 + return false;
8.471 + }
8.472 + u = _graph.source(e);
8.473 + }
8.474 + rem = (*_flow)[e];
8.475 + }
8.476 +
8.477 + _dt->addCost(n, - std::numeric_limits<Value>::max());
8.478 + _dt->addCost(n, rem);
8.479 + _dt->link(n, u);
8.480 + _dt_edges->set(n, e);
8.481 + return true;
8.482 + }
8.483 +
8.484 + void retreat(Node n) {
8.485 + _level->set(n, -1);
8.486 +
8.487 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
8.488 + Node u = _graph.target(e);
8.489 + if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
8.490 + Value rem;
8.491 + _dt->findCost(u, rem);
8.492 + _flow->set(e, rem);
8.493 + _dt->cut(u);
8.494 + _dt->addCost(u, - rem);
8.495 + _dt->addCost(u, _max_value);
8.496 + }
8.497 + }
8.498 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
8.499 + Node u = _graph.source(e);
8.500 + if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
8.501 + Value rem;
8.502 + _dt->findCost(u, rem);
8.503 + _flow->set(e, (*_capacity)[e] - rem);
8.504 + _dt->cut(u);
8.505 + _dt->addCost(u, - rem);
8.506 + _dt->addCost(u, _max_value);
8.507 + }
8.508 + }
8.509 + }
8.510 +
8.511 + void extractTrees() {
8.512 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.513 +
8.514 + Node w = _dt->findRoot(n);
8.515 +
8.516 + while (w != n) {
8.517 +
8.518 + Value rem;
8.519 + Node u = _dt->findCost(n, rem);
8.520 +
8.521 + _dt->cut(u);
8.522 + _dt->addCost(u, - rem);
8.523 + _dt->addCost(u, _max_value);
8.524 +
8.525 + Edge e = (*_dt_edges)[u];
8.526 + _dt_edges->set(u, INVALID);
8.527 +
8.528 + if (u == _graph.source(e)) {
8.529 + _flow->set(e, (*_capacity)[e] - rem);
8.530 + } else {
8.531 + _flow->set(e, rem);
8.532 + }
8.533 +
8.534 + w = _dt->findRoot(n);
8.535 + }
8.536 + }
8.537 + }
8.538 +
8.539 +
8.540 + public:
8.541 +
8.542 + /// \name Execution control The simplest way to execute the
8.543 + /// algorithm is to use the \c run() member functions.
8.544 + /// \n
8.545 + /// If you need more control on initial solution or
8.546 + /// execution then you have to call one \ref init() function and then
8.547 + /// the start() or multiple times the \c augment() member function.
8.548 +
8.549 + ///@{
8.550 +
8.551 + /// \brief Initializes the algorithm
8.552 + ///
8.553 + /// It sets the flow to empty flow.
8.554 + void init() {
8.555 + createStructures();
8.556 +
8.557 + _dt->clear();
8.558 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.559 + _dt->makeTree(n);
8.560 + _dt->addCost(n, _max_value);
8.561 + }
8.562 +
8.563 + for (EdgeIt it(_graph); it != INVALID; ++it) {
8.564 + _flow->set(it, 0);
8.565 + }
8.566 + _flow_value = 0;
8.567 + }
8.568 +
8.569 + /// \brief Initializes the algorithm
8.570 + ///
8.571 + /// Initializes the flow to the \c flowMap. The \c flowMap should
8.572 + /// contain a feasible flow, ie. in each node excluding the source
8.573 + /// and the target the incoming flow should be equal to the
8.574 + /// outgoing flow.
8.575 + template <typename FlowMap>
8.576 + void flowInit(const FlowMap& flowMap) {
8.577 + createStructures();
8.578 +
8.579 + _dt->clear();
8.580 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.581 + _dt->makeTree(n);
8.582 + _dt->addCost(n, _max_value);
8.583 + }
8.584 +
8.585 + for (EdgeIt e(_graph); e != INVALID; ++e) {
8.586 + _flow->set(e, flowMap[e]);
8.587 + }
8.588 + _flow_value = 0;
8.589 + for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
8.590 + _flow_value += (*_flow)[jt];
8.591 + }
8.592 + for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
8.593 + _flow_value -= (*_flow)[jt];
8.594 + }
8.595 + }
8.596 +
8.597 + /// \brief Initializes the algorithm
8.598 + ///
8.599 + /// Initializes the flow to the \c flowMap. The \c flowMap should
8.600 + /// contain a feasible flow, ie. in each node excluding the source
8.601 + /// and the target the incoming flow should be equal to the
8.602 + /// outgoing flow.
8.603 + /// \return %False when the given flowMap does not contain
8.604 + /// feasible flow.
8.605 + template <typename FlowMap>
8.606 + bool checkedFlowInit(const FlowMap& flowMap) {
8.607 + createStructures();
8.608 +
8.609 + _dt->clear();
8.610 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.611 + _dt->makeTree(n);
8.612 + _dt->addCost(n, _max_value);
8.613 + }
8.614 +
8.615 + for (EdgeIt e(_graph); e != INVALID; ++e) {
8.616 + _flow->set(e, flowMap[e]);
8.617 + }
8.618 + for (NodeIt it(_graph); it != INVALID; ++it) {
8.619 + if (it == _source || it == _target) continue;
8.620 + Value outFlow = 0;
8.621 + for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
8.622 + outFlow += (*_flow)[jt];
8.623 + }
8.624 + Value inFlow = 0;
8.625 + for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
8.626 + inFlow += (*_flow)[jt];
8.627 + }
8.628 + if (_tolerance.different(outFlow, inFlow)) {
8.629 + return false;
8.630 + }
8.631 + }
8.632 + for (EdgeIt it(_graph); it != INVALID; ++it) {
8.633 + if (_tolerance.less((*_flow)[it], 0)) return false;
8.634 + if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
8.635 + }
8.636 + _flow_value = 0;
8.637 + for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
8.638 + _flow_value += (*_flow)[jt];
8.639 + }
8.640 + for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
8.641 + _flow_value -= (*_flow)[jt];
8.642 + }
8.643 + return true;
8.644 + }
8.645 +
8.646 + /// \brief Executes the algorithm
8.647 + ///
8.648 + /// It runs augmenting phases by adding blocking flow until the
8.649 + /// optimal solution is reached.
8.650 + void start() {
8.651 + while (augment());
8.652 + }
8.653 +
8.654 + /// \brief Augments the flow with a blocking flow on a layered
8.655 + /// graph.
8.656 + ///
8.657 + /// This function builds a layered graph and then find a blocking
8.658 + /// flow on this graph. The number of the levels in the layered
8.659 + /// graph is strictly increasing in each augmenting phase
8.660 + /// therefore the number of the augmentings is at most \f$ n-1
8.661 + /// \f$. The length of each phase is at most \f$ O(m \log(n))
8.662 + /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
8.663 + /// \return %False when there is not residual path between the
8.664 + /// source and the target so the current flow is a feasible and
8.665 + /// optimal solution.
8.666 + bool augment() {
8.667 + Node n;
8.668 +
8.669 + if (createLayeredGraph()) {
8.670 +
8.671 + Timer bf_timer;
8.672 + initEdges();
8.673 +
8.674 + n = _dt->findRoot(_source);
8.675 + while (true) {
8.676 + Edge e;
8.677 + if (n == _target) {
8.678 + augmentPath();
8.679 + } else if (!advance(n)) {
8.680 + if (n != _source) {
8.681 + retreat(n);
8.682 + } else {
8.683 + break;
8.684 + }
8.685 + }
8.686 + n = _dt->findRoot(_source);
8.687 + }
8.688 + extractTrees();
8.689 +
8.690 + return true;
8.691 + } else {
8.692 + return false;
8.693 + }
8.694 + }
8.695 +
8.696 + /// \brief runs the algorithm.
8.697 + ///
8.698 + /// It is just a shorthand for:
8.699 + ///
8.700 + ///\code
8.701 + /// ek.init();
8.702 + /// ek.start();
8.703 + ///\endcode
8.704 + void run() {
8.705 + init();
8.706 + start();
8.707 + }
8.708 +
8.709 + /// @}
8.710 +
8.711 + /// \name Query Functions
8.712 + /// The result of the %Dijkstra algorithm can be obtained using these
8.713 + /// functions.\n
8.714 + /// Before the use of these functions,
8.715 + /// either run() or start() must be called.
8.716 +
8.717 + ///@{
8.718 +
8.719 + /// \brief Returns the value of the maximum flow.
8.720 + ///
8.721 + /// Returns the value of the maximum flow by returning the excess
8.722 + /// of the target node \c t. This value equals to the value of
8.723 + /// the maximum flow already after the first phase.
8.724 + Value flowValue() const {
8.725 + return _flow_value;
8.726 + }
8.727 +
8.728 +
8.729 + /// \brief Returns the flow on the edge.
8.730 + ///
8.731 + /// Sets the \c flowMap to the flow on the edges. This method can
8.732 + /// be called after the second phase of algorithm.
8.733 + Value flow(const Edge& edge) const {
8.734 + return (*_flow)[edge];
8.735 + }
8.736 +
8.737 + /// \brief Returns true when the node is on the source side of minimum cut.
8.738 + ///
8.739 +
8.740 + /// Returns true when the node is on the source side of minimum
8.741 + /// cut. This method can be called both after running \ref
8.742 + /// startFirstPhase() and \ref startSecondPhase().
8.743 + bool minCut(const Node& node) const {
8.744 + return (*_level)[node] == -2;
8.745 + }
8.746 +
8.747 + /// \brief Returns a minimum value cut.
8.748 + ///
8.749 + /// Sets \c cut to the characteristic vector of a minimum value cut
8.750 + /// It simply calls the minMinCut member.
8.751 + /// \retval cut Write node bool map.
8.752 + template <typename CutMap>
8.753 + void minCutMap(CutMap& cutMap) const {
8.754 + for (NodeIt n(_graph); n != INVALID; ++n) {
8.755 + cutMap.set(n, (*_level)[n] == -2);
8.756 + }
8.757 + cutMap.set(_source, true);
8.758 + }
8.759 +
8.760 + /// @}
8.761 +
8.762 + };
8.763 +}
8.764 +
8.765 +#endif
9.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
9.2 +++ b/lemon/dynamic_tree.h Sat Nov 17 20:58:11 2007 +0000
9.3 @@ -0,0 +1,1076 @@
9.4 +/* -*- C++ -*-
9.5 + *
9.6 + * This file is a part of LEMON, a generic C++ optimization library
9.7 + *
9.8 + * Copyright (C) 2003-2007
9.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
9.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
9.11 + *
9.12 + * Permission to use, modify and distribute this software is granted
9.13 + * provided that this copyright notice appears in all copies. For
9.14 + * precise terms see the accompanying LICENSE file.
9.15 + *
9.16 + * This software is provided "AS IS" with no warranty of any kind,
9.17 + * express or implied, and with no claim as to its suitability for any
9.18 + * purpose.
9.19 + *
9.20 + */
9.21 +#ifndef LEMON_DYNAMIC_TREE_H
9.22 +#define LEMON_DYNAMIC_TREE_H
9.23 +
9.24 +/// \ingroup auxdata
9.25 +/// \file
9.26 +/// \brief The dynamic tree data structure of Sleator and Tarjan.
9.27 +///
9.28 +
9.29 +#include <vector>
9.30 +#include <limits>
9.31 +#include <lemon/tolerance.h>
9.32 +
9.33 +namespace lemon {
9.34 +
9.35 + /// \ingroup auxdata
9.36 + ///
9.37 + /// \brief The dynamic tree data structure of Sleator and Tarjan.
9.38 + ///
9.39 + /// This class provides an implementation of the dynamic tree data
9.40 + /// structure for maintaining a set of node-disjoint rooted
9.41 + /// trees. Each item has an associated value, and the item with
9.42 + /// minimum value can be find in \f$O(\log(n)\f$ on the path from a
9.43 + /// node to the its root, and the items on such path can be
9.44 + /// increased or decreased globally with a certain value in the same
9.45 + /// running time. We regard a tree edge as directed toward the root,
9.46 + /// that is from child to parent. Its structure can be modified by
9.47 + /// two basic operations: \e link(v,w) adds an edge between a root v
9.48 + /// and a node w in a different component; \e cut(v) removes the
9.49 + /// edge between v and its parent.
9.50 + ///
9.51 + /// \param _Value The value type of the items.
9.52 + /// \param _ItemIntMap Converts item type of node to integer.
9.53 + /// \param _Tolerance The tolerance class to handle computation
9.54 + /// problems.
9.55 + /// \param _enableSize If true then the data structre manatain the
9.56 + /// size of each tree. The feature is used in \ref GoldbergTarjan
9.57 + /// algorithm. The default value is true.
9.58 + ///
9.59 + /// \author Hamori Tamas
9.60 +#ifdef DOXYGEN
9.61 + template <typename _Value, typename _ItemIntMap,
9.62 + typename _Tolerance, bool _enableSize>
9.63 +#else
9.64 + template <typename _Value, typename _ItemIntMap,
9.65 + typename _Tolerance = lemon::Tolerance<_Value>,
9.66 + bool _enableSize = true>
9.67 +#endif
9.68 + class DynamicTree {
9.69 + public:
9.70 + /// \brief The integer map on the items.
9.71 + typedef _ItemIntMap ItemIntMap;
9.72 + /// \brief The item type of nodes.
9.73 + typedef typename ItemIntMap::Key Item;
9.74 + /// \brief The value type of the algorithms.
9.75 + typedef _Value Value;
9.76 + /// \brief The tolerance used by the algorithm.
9.77 + typedef _Tolerance Tolerance;
9.78 +
9.79 + private:
9.80 +
9.81 + class ItemData;
9.82 +
9.83 + std::vector<ItemData> _data;
9.84 + ItemIntMap &_iim;
9.85 + Value _max_value;
9.86 + Tolerance _tolerance;
9.87 +
9.88 + public:
9.89 + /// \brief The constructor of the class.
9.90 + ///
9.91 + /// \param iim The integer map on the items.
9.92 + /// \param tolerance Tolerance class.
9.93 + DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
9.94 + : _iim(iim), _max_value(std::numeric_limits<Value>::max()),
9.95 + _tolerance(tolerance) {}
9.96 +
9.97 + ~DynamicTree() {}
9.98 +
9.99 + /// \brief Clears the data structure
9.100 + ///
9.101 + /// Clears the data structure
9.102 + void clear() {
9.103 + _data.clear();
9.104 + }
9.105 +
9.106 + /// \brief Sets the tolerance used by algorithm.
9.107 + ///
9.108 + /// Sets the tolerance used by algorithm.
9.109 + void tolerance(const Tolerance& tolerance) const {
9.110 + _tolerance = tolerance;
9.111 + return *this;
9.112 + }
9.113 +
9.114 + /// \brief Returns the tolerance used by algorithm.
9.115 + ///
9.116 + /// Returns the tolerance used by algorithm.
9.117 + const Tolerance& tolerance() const {
9.118 + return tolerance;
9.119 + }
9.120 +
9.121 + /// \brief Create a new tree containing a single node with cost zero.
9.122 + void makeTree(const Item &item) {
9.123 + _data[makePath(item)].successor = -1;
9.124 + }
9.125 +
9.126 + /// \brief Return the root of the tree containing node with itemtype
9.127 + /// \e item.
9.128 + Item findRoot(const Item &item) {
9.129 + return _data[findTail(expose(_iim[item]))].id;
9.130 + }
9.131 +
9.132 + /// \brief Return the the value of nodes in the tree containing
9.133 + /// node with itemtype \e item.
9.134 + int findSize(const Item &item) {
9.135 + return _data[expose(_iim[item])].size;
9.136 + }
9.137 +
9.138 + /// \brief Return the minimum cost containing node.
9.139 + ///
9.140 + /// Return into \e d the minimum cost on the tree path from \e item
9.141 + /// to findRoot(item). Return the last item (closest to its root)
9.142 + /// on this path of the minimum cost.
9.143 + Item findCost(const Item &item, Value& d){
9.144 + return _data[findPathCost(expose(_iim[item]),d)].id;
9.145 + }
9.146 +
9.147 + /// \brief Add \e x value to the cost of every node on the path from
9.148 + /// \e item to findRoot(item).
9.149 + void addCost(const Item &item, Value x){
9.150 + addPathCost(expose(_iim[item]), x);
9.151 + }
9.152 +
9.153 + /// \brief Combine the trees containing nodes \e item1 and \e item2
9.154 + /// by adding an edge from \e item1 \e item2.
9.155 + ///
9.156 + /// This operation assumes that \e item1 is root and \e item2 is in
9.157 + /// a different tree.
9.158 + void link(const Item &item1, const Item &item2){
9.159 + int v = _iim[item1];
9.160 + int w = _iim[item2];
9.161 + int p = expose(w);
9.162 + join(-1, expose(v), p);
9.163 + _data[v].successor = -1;
9.164 + _data[v].size += _data[p].size;
9.165 +
9.166 + }
9.167 +
9.168 + /// \brief Divide the tree containing node \e item into two trees by
9.169 + /// deleting the edge out of \e item.
9.170 + ///
9.171 + /// This operation assumes that \e item is not a tree root.
9.172 + void cut(const Item &item) {
9.173 + int v = _iim[item];
9.174 + int p, q;
9.175 + expose(v);
9.176 + split(p, v, q);
9.177 + if (p != -1) {
9.178 + _data[p].successor = v;
9.179 + }
9.180 + _data[v].size -= _data[q].size;
9.181 + if (q != -1) {
9.182 + _data[q].successor = _data[v].successor;
9.183 + }
9.184 + _data[v].successor = -1;
9.185 + }
9.186 +
9.187 + ///\brief
9.188 + Item parent(const Item &item){
9.189 + return _data[_iim[item].p].id;
9.190 + }
9.191 +
9.192 + ///\brief Return the upper bound of the costs.
9.193 + Value maxValue() const {
9.194 + return _max_value;
9.195 + }
9.196 +
9.197 + private:
9.198 +
9.199 + int makePath(const Item &item) {
9.200 + _iim.set(item, _data.size());
9.201 + ItemData v(item);
9.202 + _data.push_back(v);
9.203 + return _iim[item];
9.204 + }
9.205 +
9.206 + int findPath(int v){
9.207 + splay(v);
9.208 + return v;
9.209 + }
9.210 +
9.211 + int findPathCost(int p, Value &d){
9.212 + while ((_data[p].right != -1 &&
9.213 + !_tolerance.less(0, _data[_data[p].right].dmin)) ||
9.214 + (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
9.215 + if (_data[p].right != -1 &&
9.216 + !_tolerance.less(0, _data[_data[p].right].dmin)) {
9.217 + p = _data[p].right;
9.218 + } else if (_data[p].left != -1 &&
9.219 + !_tolerance.less(0, _data[_data[p].left].dmin)){
9.220 + p = _data[p].left;
9.221 + }
9.222 + }
9.223 + splay(p);
9.224 + d = _data[p].dmin;
9.225 + return p;
9.226 + }
9.227 +
9.228 + int findTail(int p){
9.229 + while (_data[p].right != -1) {
9.230 + p = _data[p].right;
9.231 + }
9.232 + splay(p);
9.233 + return p;
9.234 + }
9.235 +
9.236 + void addPathCost(int p, Value x){
9.237 + if (!_tolerance.less(x, _max_value)) {
9.238 + _data[p].dmin = x;_data[p].dcost = x;
9.239 + } else if (!_tolerance.less(-x, _max_value)) {
9.240 + _data[p].dmin = 0;
9.241 + _data[p].dcost = 0;
9.242 + } else {
9.243 + _data[p].dmin += x;
9.244 + }
9.245 + }
9.246 +
9.247 + void join(int p, int v, int q) {
9.248 + Value min = _max_value;
9.249 + Value pmin = _max_value;
9.250 + Value vmin = _data[v].dmin;
9.251 + Value qmin = _max_value;
9.252 + if (p != -1){
9.253 + pmin = _data[p].dmin;
9.254 + }
9.255 + if (q != -1){
9.256 + qmin = _data[q].dmin;
9.257 + }
9.258 +
9.259 + if (_tolerance.less(vmin, qmin)) {
9.260 + if (_tolerance.less(vmin,pmin)) {
9.261 + min = vmin;
9.262 + } else {
9.263 + min = pmin;
9.264 + }
9.265 + } else if (_tolerance.less(qmin,pmin)) {
9.266 + min = qmin;
9.267 + } else {
9.268 + min = pmin;
9.269 + }
9.270 +
9.271 + if (p != -1){
9.272 + _data[p].parent = v;
9.273 + _data[p].dmin -= min;
9.274 + }
9.275 + if (q!=-1){
9.276 + _data[q].parent = v;
9.277 + if (_tolerance.less(_data[q].dmin,_max_value)) {
9.278 + _data[q].dmin -= min;
9.279 + }
9.280 + }
9.281 + _data[v].left = p;
9.282 + _data[v].right = q;
9.283 + if (_tolerance.less(min,_max_value)) {
9.284 + _data[v].dcost = _data[v].dmin - min;
9.285 + }
9.286 + _data[v].dmin = min;
9.287 + }
9.288 +
9.289 + void split(int &p, int v, int &q){
9.290 + splay(v);
9.291 + p = -1;
9.292 + if (_data[v].left != -1){
9.293 + p = _data[v].left;
9.294 + _data[p].dmin += _data[v].dmin;
9.295 + _data[p].parent = -1;
9.296 + _data[v].left = -1;
9.297 + }
9.298 + q = -1;
9.299 + if (_data[v].right != -1) {
9.300 + q=_data[v].right;
9.301 + if (_tolerance.less(_data[q].dmin, _max_value)) {
9.302 + _data[q].dmin += _data[v].dmin;
9.303 + }
9.304 + _data[q].parent = -1;
9.305 + _data[v].right = -1;
9.306 + }
9.307 + if (_tolerance.less(_data[v].dcost, _max_value)) {
9.308 + _data[v].dmin += _data[v].dcost;
9.309 + _data[v].dcost = 0;
9.310 + } else {
9.311 + _data[v].dmin = _data[v].dcost;
9.312 + }
9.313 + }
9.314 +
9.315 + int expose(int v) {
9.316 + int p, q, r, w;
9.317 + p = -1;
9.318 + while (v != -1) {
9.319 + w = _data[findPath(v)].successor;
9.320 + split(q, v, r);
9.321 + if (q != -1) {
9.322 + _data[q].successor = v;
9.323 + }
9.324 + join(p, v, r);
9.325 + p = v;
9.326 + v = w;
9.327 + }
9.328 + _data[p].successor = -1;
9.329 + return p;
9.330 + }
9.331 +
9.332 + void splay(int v) {
9.333 + while (_data[v].parent != -1) {
9.334 + if (v == _data[_data[v].parent].left) {
9.335 + if (_data[_data[v].parent].parent == -1) {
9.336 + zig(v);
9.337 + } else {
9.338 + if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
9.339 + zig(_data[v].parent);
9.340 + zig(v);
9.341 + } else {
9.342 + zig(v);
9.343 + zag(v);
9.344 + }
9.345 + }
9.346 + } else {
9.347 + if (_data[_data[v].parent].parent == -1) {
9.348 + zag(v);
9.349 + } else {
9.350 + if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
9.351 + zag(v);
9.352 + zig(v);
9.353 + } else {
9.354 + zag(_data[v].parent);
9.355 + zag(v);
9.356 + }
9.357 + }
9.358 + }
9.359 + }
9.360 + }
9.361 +
9.362 +
9.363 + void zig(int v) {
9.364 + Value min = _data[_data[v].parent].dmin;
9.365 + int a = _data[v].parent;
9.366 +
9.367 + Value aa = _data[a].dcost;
9.368 + if (_tolerance.less(aa, _max_value)) {
9.369 + aa+= min;
9.370 + }
9.371 +
9.372 +
9.373 + int b = v;
9.374 + Value ab = min + _data[b].dmin;
9.375 + Value bb = _data[b].dcost;
9.376 + if (_tolerance.less(bb, _max_value)) {
9.377 + bb+= ab;
9.378 + }
9.379 +
9.380 + int c = -1;
9.381 + Value cc = _max_value;
9.382 + if (_data[a].right != -1) {
9.383 + c = _data[a].right;
9.384 + cc = _data[c].dmin;
9.385 + if (_tolerance.less(cc, _max_value)) {
9.386 + cc+=min;
9.387 + }
9.388 + }
9.389 +
9.390 + int d = -1;
9.391 + Value dd = _max_value;
9.392 + if (_data[v].left != -1){
9.393 + d = _data[v].left;
9.394 + dd = ab + _data[d].dmin;
9.395 + }
9.396 +
9.397 + int e = -1;
9.398 + Value ee = _max_value;
9.399 + if (_data[v].right != -1) {
9.400 + e = _data[v].right;
9.401 + ee = ab + _data[e].dmin;
9.402 + }
9.403 +
9.404 + Value min2;
9.405 + if (_tolerance.less(0, _data[b].dmin) ||
9.406 + (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
9.407 + min2 = min;
9.408 + } else {
9.409 + if (_tolerance.less(aa, cc)) {
9.410 + if (_tolerance.less(aa, ee)) {
9.411 + min2 = aa;
9.412 + } else {
9.413 + min2 = ee;
9.414 + }
9.415 + } else if (_tolerance.less(cc, ee)) {
9.416 + min2 = cc;
9.417 + } else {
9.418 + min2 = ee;
9.419 + }
9.420 + }
9.421 +
9.422 + _data[a].dcost = aa;
9.423 + if (_tolerance.less(aa, _max_value)) {
9.424 + _data[a].dcost -= min2;
9.425 + }
9.426 + _data[a].dmin = min2;
9.427 + if (_tolerance.less(min2,_max_value)) {
9.428 + _data[a].dmin -= min;
9.429 + }
9.430 + _data[a].size -= _data[b].size;
9.431 + _data[b].dcost = bb;
9.432 + if (_tolerance.less(bb, _max_value)) {
9.433 + _data[b].dcost -= min;
9.434 + }
9.435 + _data[b].dmin = min;
9.436 + _data[b].size += _data[a].size;
9.437 + if (c != -1) {
9.438 + _data[c].dmin = cc;
9.439 + if (_tolerance.less(cc, _max_value)) {
9.440 + _data[c].dmin -= min2;
9.441 + }
9.442 + }
9.443 + if (d != -1) {
9.444 + _data[d].dmin = dd - min;
9.445 + _data[a].size += _data[d].size;
9.446 + _data[b].size -= _data[d].size;
9.447 + }
9.448 + if (e != -1) {
9.449 + _data[e].dmin = ee - min2;
9.450 + }
9.451 +
9.452 + int w = _data[v].parent;
9.453 + _data[v].successor = _data[w].successor;
9.454 + _data[w].successor = -1;
9.455 + _data[v].parent = _data[w].parent;
9.456 + _data[w].parent = v;
9.457 + _data[w].left = _data[v].right;
9.458 + _data[v].right = w;
9.459 + if (_data[v].parent != -1){
9.460 + if (_data[_data[v].parent].right == w) {
9.461 + _data[_data[v].parent].right = v;
9.462 + } else {
9.463 + _data[_data[v].parent].left = v;
9.464 + }
9.465 + }
9.466 + if (_data[w].left != -1){
9.467 + _data[_data[w].left].parent = w;
9.468 + }
9.469 + }
9.470 +
9.471 +
9.472 + void zag(int v) {
9.473 +
9.474 + Value min = _data[_data[v].parent].dmin;
9.475 +
9.476 + int a = _data[v].parent;
9.477 + Value aa = _data[a].dcost;
9.478 + if (_tolerance.less(aa, _max_value)) {
9.479 + aa += min;
9.480 + }
9.481 +
9.482 + int b = v;
9.483 + Value ab = min + _data[b].dmin;
9.484 + Value bb = _data[b].dcost;
9.485 + if (_tolerance.less(bb, _max_value)) {
9.486 + bb += ab;
9.487 + }
9.488 +
9.489 + int c = -1;
9.490 + Value cc = _max_value;
9.491 + if (_data[a].left != -1){
9.492 + c = _data[a].left;
9.493 + cc = min + _data[c].dmin;
9.494 + }
9.495 +
9.496 + int d = -1;
9.497 + Value dd = _max_value;
9.498 + if (_data[v].right!=-1) {
9.499 + d = _data[v].right;
9.500 + dd = _data[d].dmin;
9.501 + if (_tolerance.less(dd, _max_value)) {
9.502 + dd += ab;
9.503 + }
9.504 + }
9.505 +
9.506 + int e = -1;
9.507 + Value ee = _max_value;
9.508 + if (_data[v].left != -1){
9.509 + e = _data[v].left;
9.510 + ee = ab + _data[e].dmin;
9.511 + }
9.512 +
9.513 + Value min2;
9.514 + if (_tolerance.less(0, _data[b].dmin) ||
9.515 + (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
9.516 + min2 = min;
9.517 + } else {
9.518 + if (_tolerance.less(aa, cc)) {
9.519 + if (_tolerance.less(aa, ee)) {
9.520 + min2 = aa;
9.521 + } else {
9.522 + min2 = ee;
9.523 + }
9.524 + } else if (_tolerance.less(cc, ee)) {
9.525 + min2 = cc;
9.526 + } else {
9.527 + min2 = ee;
9.528 + }
9.529 + }
9.530 + _data[a].dcost = aa;
9.531 + if (_tolerance.less(aa, _max_value)) {
9.532 + _data[a].dcost -= min2;
9.533 + }
9.534 + _data[a].dmin = min2;
9.535 + if (_tolerance.less(min2, _max_value)) {
9.536 + _data[a].dmin -= min;
9.537 + }
9.538 + _data[a].size -= _data[b].size;
9.539 + _data[b].dcost = bb;
9.540 + if (_tolerance.less(bb, _max_value)) {
9.541 + _data[b].dcost -= min;
9.542 + }
9.543 + _data[b].dmin = min;
9.544 + _data[b].size += _data[a].size;
9.545 + if (c != -1) {
9.546 + _data[c].dmin = cc - min2;
9.547 + }
9.548 + if (d != -1) {
9.549 + _data[d].dmin = dd;
9.550 + _data[a].size += _data[d].size;
9.551 + _data[b].size -= _data[d].size;
9.552 + if (_tolerance.less(dd, _max_value)) {
9.553 + _data[d].dmin -= min;
9.554 + }
9.555 + }
9.556 + if (e != -1) {
9.557 + _data[e].dmin = ee - min2;
9.558 + }
9.559 +
9.560 + int w = _data[v].parent;
9.561 + _data[v].successor = _data[w].successor;
9.562 + _data[w].successor = -1;
9.563 + _data[v].parent = _data[w].parent;
9.564 + _data[w].parent = v;
9.565 + _data[w].right = _data[v].left;
9.566 + _data[v].left = w;
9.567 + if (_data[v].parent != -1){
9.568 + if (_data[_data[v].parent].left == w) {
9.569 + _data[_data[v].parent].left = v;
9.570 + } else {
9.571 + _data[_data[v].parent].right = v;
9.572 + }
9.573 + }
9.574 + if (_data[w].right != -1){
9.575 + _data[_data[w].right].parent = w;
9.576 + }
9.577 + }
9.578 +
9.579 + private:
9.580 +
9.581 + class ItemData {
9.582 + public:
9.583 + Item id;
9.584 + int size;
9.585 + int successor;
9.586 + int parent;
9.587 + int left;
9.588 + int right;
9.589 + Value dmin;
9.590 + Value dcost;
9.591 +
9.592 + public:
9.593 + ItemData(const Item &item)
9.594 + : id(item), size(1), successor(), parent(-1),
9.595 + left(-1), right(-1), dmin(0), dcost(0) {}
9.596 + };
9.597 +
9.598 + };
9.599 +
9.600 + template <typename _Value, typename _ItemIntMap, typename _Tolerance>
9.601 + class DynamicTree<_Value, _ItemIntMap, _Tolerance, false> {
9.602 + public:
9.603 + typedef _ItemIntMap ItemIntMap;
9.604 + typedef typename ItemIntMap::Key Item;
9.605 + typedef _Value Value;
9.606 + typedef _Tolerance Tolerance;
9.607 +
9.608 + private:
9.609 +
9.610 + class ItemData;
9.611 +
9.612 + std::vector<ItemData> _data;
9.613 + ItemIntMap &_iim;
9.614 + Value _max_value;
9.615 + Tolerance _tolerance;
9.616 +
9.617 + public:
9.618 + DynamicTree(ItemIntMap &iim, const Tolerance& tolerance = Tolerance())
9.619 + : _iim(iim), _max_value(std::numeric_limits<Value>::max()),
9.620 + _tolerance(tolerance) {}
9.621 +
9.622 + ~DynamicTree() {}
9.623 +
9.624 + void clear() {
9.625 + _data.clear();
9.626 + }
9.627 +
9.628 + void tolerance(const Tolerance& tolerance) const {
9.629 + _tolerance = tolerance;
9.630 + return *this;
9.631 + }
9.632 +
9.633 + const Tolerance& tolerance() const {
9.634 + return tolerance;
9.635 + }
9.636 +
9.637 + void makeTree(const Item &item) {
9.638 + _data[makePath(item)].successor = -1;
9.639 + }
9.640 +
9.641 + Item findRoot(const Item &item) {
9.642 + return _data[findTail(expose(_iim[item]))].id;
9.643 + }
9.644 +
9.645 + Item findCost(const Item &item, Value& d){
9.646 + return _data[findPathCost(expose(_iim[item]),d)].id;
9.647 + }
9.648 +
9.649 + void addCost(const Item &item, Value x){
9.650 + addPathCost(expose(_iim[item]), x);
9.651 + }
9.652 +
9.653 + void link(const Item &item1, const Item &item2){
9.654 + int v = _iim[item1];
9.655 + int w = _iim[item2];
9.656 + int p = expose(w);
9.657 + join(-1, expose(v), p);
9.658 + _data[v].successor = -1;
9.659 + }
9.660 +
9.661 + void cut(const Item &item) {
9.662 + int v = _iim[item];
9.663 + int p, q;
9.664 + expose(v);
9.665 + split(p, v, q);
9.666 + if (p != -1) {
9.667 + _data[p].successor = v;
9.668 + }
9.669 + if (q != -1) {
9.670 + _data[q].successor = _data[v].successor;
9.671 + }
9.672 + _data[v].successor = -1;
9.673 + }
9.674 +
9.675 + Item parent(const Item &item){
9.676 + return _data[_iim[item].p].id;
9.677 + }
9.678 +
9.679 + Value maxValue() const {
9.680 + return _max_value;
9.681 + }
9.682 +
9.683 + private:
9.684 +
9.685 + int makePath(const Item &item) {
9.686 + _iim.set(item, _data.size());
9.687 + ItemData v(item);
9.688 + _data.push_back(v);
9.689 + return _iim[item];
9.690 + }
9.691 +
9.692 + int findPath(int v){
9.693 + splay(v);
9.694 + return v;
9.695 + }
9.696 +
9.697 + int findPathCost(int p, Value &d){
9.698 + while ((_data[p].right != -1 &&
9.699 + !_tolerance.less(0, _data[_data[p].right].dmin)) ||
9.700 + (_data[p].left != -1 && _tolerance.less(0, _data[p].dcost))) {
9.701 + if (_data[p].right != -1 &&
9.702 + !_tolerance.less(0, _data[_data[p].right].dmin)) {
9.703 + p = _data[p].right;
9.704 + } else if (_data[p].left != -1 &&
9.705 + !_tolerance.less(0, _data[_data[p].left].dmin)){
9.706 + p = _data[p].left;
9.707 + }
9.708 + }
9.709 + splay(p);
9.710 + d = _data[p].dmin;
9.711 + return p;
9.712 + }
9.713 +
9.714 + int findTail(int p){
9.715 + while (_data[p].right != -1) {
9.716 + p = _data[p].right;
9.717 + }
9.718 + splay(p);
9.719 + return p;
9.720 + }
9.721 +
9.722 + void addPathCost(int p, Value x){
9.723 + if (!_tolerance.less(x, _max_value)) {
9.724 + _data[p].dmin = x;_data[p].dcost = x;
9.725 + } else if (!_tolerance.less(-x, _max_value)) {
9.726 + _data[p].dmin = 0;
9.727 + _data[p].dcost = 0;
9.728 + } else {
9.729 + _data[p].dmin += x;
9.730 + }
9.731 + }
9.732 +
9.733 + void join(int p, int v, int q) {
9.734 + Value min = _max_value;
9.735 + Value pmin = _max_value;
9.736 + Value vmin = _data[v].dmin;
9.737 + Value qmin = _max_value;
9.738 + if (p != -1){
9.739 + pmin = _data[p].dmin;
9.740 + }
9.741 + if (q != -1){
9.742 + qmin = _data[q].dmin;
9.743 + }
9.744 +
9.745 + if (_tolerance.less(vmin, qmin)) {
9.746 + if (_tolerance.less(vmin,pmin)) {
9.747 + min = vmin;
9.748 + } else {
9.749 + min = pmin;
9.750 + }
9.751 + } else if (_tolerance.less(qmin,pmin)) {
9.752 + min = qmin;
9.753 + } else {
9.754 + min = pmin;
9.755 + }
9.756 +
9.757 + if (p != -1){
9.758 + _data[p].parent = v;
9.759 + _data[p].dmin -= min;
9.760 + }
9.761 + if (q!=-1){
9.762 + _data[q].parent = v;
9.763 + if (_tolerance.less(_data[q].dmin,_max_value)) {
9.764 + _data[q].dmin -= min;
9.765 + }
9.766 + }
9.767 + _data[v].left = p;
9.768 + _data[v].right = q;
9.769 + if (_tolerance.less(min,_max_value)) {
9.770 + _data[v].dcost = _data[v].dmin - min;
9.771 + }
9.772 + _data[v].dmin = min;
9.773 + }
9.774 +
9.775 + void split(int &p, int v, int &q){
9.776 + splay(v);
9.777 + p = -1;
9.778 + if (_data[v].left != -1){
9.779 + p = _data[v].left;
9.780 + _data[p].dmin += _data[v].dmin;
9.781 + _data[p].parent = -1;
9.782 + _data[v].left = -1;
9.783 + }
9.784 + q = -1;
9.785 + if (_data[v].right != -1) {
9.786 + q=_data[v].right;
9.787 + if (_tolerance.less(_data[q].dmin, _max_value)) {
9.788 + _data[q].dmin += _data[v].dmin;
9.789 + }
9.790 + _data[q].parent = -1;
9.791 + _data[v].right = -1;
9.792 + }
9.793 + if (_tolerance.less(_data[v].dcost, _max_value)) {
9.794 + _data[v].dmin += _data[v].dcost;
9.795 + _data[v].dcost = 0;
9.796 + } else {
9.797 + _data[v].dmin = _data[v].dcost;
9.798 + }
9.799 + }
9.800 +
9.801 + int expose(int v) {
9.802 + int p, q, r, w;
9.803 + p = -1;
9.804 + while (v != -1) {
9.805 + w = _data[findPath(v)].successor;
9.806 + split(q, v, r);
9.807 + if (q != -1) {
9.808 + _data[q].successor = v;
9.809 + }
9.810 + join(p, v, r);
9.811 + p = v;
9.812 + v = w;
9.813 + }
9.814 + _data[p].successor = -1;
9.815 + return p;
9.816 + }
9.817 +
9.818 + void splay(int v) {
9.819 + while (_data[v].parent != -1) {
9.820 + if (v == _data[_data[v].parent].left) {
9.821 + if (_data[_data[v].parent].parent == -1) {
9.822 + zig(v);
9.823 + } else {
9.824 + if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
9.825 + zig(_data[v].parent);
9.826 + zig(v);
9.827 + } else {
9.828 + zig(v);
9.829 + zag(v);
9.830 + }
9.831 + }
9.832 + } else {
9.833 + if (_data[_data[v].parent].parent == -1) {
9.834 + zag(v);
9.835 + } else {
9.836 + if (_data[v].parent == _data[_data[_data[v].parent].parent].left) {
9.837 + zag(v);
9.838 + zig(v);
9.839 + } else {
9.840 + zag(_data[v].parent);
9.841 + zag(v);
9.842 + }
9.843 + }
9.844 + }
9.845 + }
9.846 + }
9.847 +
9.848 +
9.849 + void zig(int v) {
9.850 + Value min = _data[_data[v].parent].dmin;
9.851 + int a = _data[v].parent;
9.852 +
9.853 + Value aa = _data[a].dcost;
9.854 + if (_tolerance.less(aa, _max_value)) {
9.855 + aa+= min;
9.856 + }
9.857 +
9.858 +
9.859 + int b = v;
9.860 + Value ab = min + _data[b].dmin;
9.861 + Value bb = _data[b].dcost;
9.862 + if (_tolerance.less(bb, _max_value)) {
9.863 + bb+= ab;
9.864 + }
9.865 +
9.866 + int c = -1;
9.867 + Value cc = _max_value;
9.868 + if (_data[a].right != -1) {
9.869 + c = _data[a].right;
9.870 + cc = _data[c].dmin;
9.871 + if (_tolerance.less(cc, _max_value)) {
9.872 + cc+=min;
9.873 + }
9.874 + }
9.875 +
9.876 + int d = -1;
9.877 + Value dd = _max_value;
9.878 + if (_data[v].left != -1){
9.879 + d = _data[v].left;
9.880 + dd = ab + _data[d].dmin;
9.881 + }
9.882 +
9.883 + int e = -1;
9.884 + Value ee = _max_value;
9.885 + if (_data[v].right != -1) {
9.886 + e = _data[v].right;
9.887 + ee = ab + _data[e].dmin;
9.888 + }
9.889 +
9.890 + Value min2;
9.891 + if (_tolerance.less(0, _data[b].dmin) ||
9.892 + (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
9.893 + min2 = min;
9.894 + } else {
9.895 + if (_tolerance.less(aa, cc)) {
9.896 + if (_tolerance.less(aa, ee)) {
9.897 + min2 = aa;
9.898 + } else {
9.899 + min2 = ee;
9.900 + }
9.901 + } else if (_tolerance.less(cc, ee)) {
9.902 + min2 = cc;
9.903 + } else {
9.904 + min2 = ee;
9.905 + }
9.906 + }
9.907 +
9.908 + _data[a].dcost = aa;
9.909 + if (_tolerance.less(aa, _max_value)) {
9.910 + _data[a].dcost -= min2;
9.911 + }
9.912 + _data[a].dmin = min2;
9.913 + if (_tolerance.less(min2,_max_value)) {
9.914 + _data[a].dmin -= min;
9.915 + }
9.916 + _data[b].dcost = bb;
9.917 + if (_tolerance.less(bb, _max_value)) {
9.918 + _data[b].dcost -= min;
9.919 + }
9.920 + _data[b].dmin = min;
9.921 + if (c != -1) {
9.922 + _data[c].dmin = cc;
9.923 + if (_tolerance.less(cc, _max_value)) {
9.924 + _data[c].dmin -= min2;
9.925 + }
9.926 + }
9.927 + if (d != -1) {
9.928 + _data[d].dmin = dd - min;
9.929 + }
9.930 + if (e != -1) {
9.931 + _data[e].dmin = ee - min2;
9.932 + }
9.933 +
9.934 + int w = _data[v].parent;
9.935 + _data[v].successor = _data[w].successor;
9.936 + _data[w].successor = -1;
9.937 + _data[v].parent = _data[w].parent;
9.938 + _data[w].parent = v;
9.939 + _data[w].left = _data[v].right;
9.940 + _data[v].right = w;
9.941 + if (_data[v].parent != -1){
9.942 + if (_data[_data[v].parent].right == w) {
9.943 + _data[_data[v].parent].right = v;
9.944 + } else {
9.945 + _data[_data[v].parent].left = v;
9.946 + }
9.947 + }
9.948 + if (_data[w].left != -1){
9.949 + _data[_data[w].left].parent = w;
9.950 + }
9.951 + }
9.952 +
9.953 +
9.954 + void zag(int v) {
9.955 +
9.956 + Value min = _data[_data[v].parent].dmin;
9.957 +
9.958 + int a = _data[v].parent;
9.959 + Value aa = _data[a].dcost;
9.960 + if (_tolerance.less(aa, _max_value)) {
9.961 + aa += min;
9.962 + }
9.963 +
9.964 + int b = v;
9.965 + Value ab = min + _data[b].dmin;
9.966 + Value bb = _data[b].dcost;
9.967 + if (_tolerance.less(bb, _max_value)) {
9.968 + bb += ab;
9.969 + }
9.970 +
9.971 + int c = -1;
9.972 + Value cc = _max_value;
9.973 + if (_data[a].left != -1){
9.974 + c = _data[a].left;
9.975 + cc = min + _data[c].dmin;
9.976 + }
9.977 +
9.978 + int d = -1;
9.979 + Value dd = _max_value;
9.980 + if (_data[v].right!=-1) {
9.981 + d = _data[v].right;
9.982 + dd = _data[d].dmin;
9.983 + if (_tolerance.less(dd, _max_value)) {
9.984 + dd += ab;
9.985 + }
9.986 + }
9.987 +
9.988 + int e = -1;
9.989 + Value ee = _max_value;
9.990 + if (_data[v].left != -1){
9.991 + e = _data[v].left;
9.992 + ee = ab + _data[e].dmin;
9.993 + }
9.994 +
9.995 + Value min2;
9.996 + if (_tolerance.less(0, _data[b].dmin) ||
9.997 + (e != -1 && !_tolerance.less(0, _data[e].dmin))) {
9.998 + min2 = min;
9.999 + } else {
9.1000 + if (_tolerance.less(aa, cc)) {
9.1001 + if (_tolerance.less(aa, ee)) {
9.1002 + min2 = aa;
9.1003 + } else {
9.1004 + min2 = ee;
9.1005 + }
9.1006 + } else if (_tolerance.less(cc, ee)) {
9.1007 + min2 = cc;
9.1008 + } else {
9.1009 + min2 = ee;
9.1010 + }
9.1011 + }
9.1012 + _data[a].dcost = aa;
9.1013 + if (_tolerance.less(aa, _max_value)) {
9.1014 + _data[a].dcost -= min2;
9.1015 + }
9.1016 + _data[a].dmin = min2;
9.1017 + if (_tolerance.less(min2, _max_value)) {
9.1018 + _data[a].dmin -= min;
9.1019 + }
9.1020 + _data[b].dcost = bb;
9.1021 + if (_tolerance.less(bb, _max_value)) {
9.1022 + _data[b].dcost -= min;
9.1023 + }
9.1024 + _data[b].dmin = min;
9.1025 + if (c != -1) {
9.1026 + _data[c].dmin = cc - min2;
9.1027 + }
9.1028 + if (d != -1) {
9.1029 + _data[d].dmin = dd;
9.1030 + if (_tolerance.less(dd, _max_value)) {
9.1031 + _data[d].dmin -= min;
9.1032 + }
9.1033 + }
9.1034 + if (e != -1) {
9.1035 + _data[e].dmin = ee - min2;
9.1036 + }
9.1037 +
9.1038 + int w = _data[v].parent;
9.1039 + _data[v].successor = _data[w].successor;
9.1040 + _data[w].successor = -1;
9.1041 + _data[v].parent = _data[w].parent;
9.1042 + _data[w].parent = v;
9.1043 + _data[w].right = _data[v].left;
9.1044 + _data[v].left = w;
9.1045 + if (_data[v].parent != -1){
9.1046 + if (_data[_data[v].parent].left == w) {
9.1047 + _data[_data[v].parent].left = v;
9.1048 + } else {
9.1049 + _data[_data[v].parent].right = v;
9.1050 + }
9.1051 + }
9.1052 + if (_data[w].right != -1){
9.1053 + _data[_data[w].right].parent = w;
9.1054 + }
9.1055 + }
9.1056 +
9.1057 + private:
9.1058 +
9.1059 + class ItemData {
9.1060 + public:
9.1061 + Item id;
9.1062 + int successor;
9.1063 + int parent;
9.1064 + int left;
9.1065 + int right;
9.1066 + Value dmin;
9.1067 + Value dcost;
9.1068 +
9.1069 + public:
9.1070 + ItemData(const Item &item)
9.1071 + : id(item), successor(), parent(-1),
9.1072 + left(-1), right(-1), dmin(0), dcost(0) {}
9.1073 + };
9.1074 +
9.1075 + };
9.1076 +
9.1077 +}
9.1078 +
9.1079 +#endif
10.1 --- a/lemon/edmonds_karp.h Wed Nov 14 17:53:08 2007 +0000
10.2 +++ b/lemon/edmonds_karp.h Sat Nov 17 20:58:11 2007 +0000
10.3 @@ -23,47 +23,95 @@
10.4 /// \ingroup max_flow
10.5 /// \brief Implementation of the Edmonds-Karp algorithm.
10.6
10.7 -#include <lemon/graph_adaptor.h>
10.8 #include <lemon/tolerance.h>
10.9 -#include <lemon/bfs.h>
10.10 +#include <vector>
10.11
10.12 namespace lemon {
10.13
10.14 + /// \brief Default traits class of EdmondsKarp class.
10.15 + ///
10.16 + /// Default traits class of EdmondsKarp class.
10.17 + /// \param _Graph Graph type.
10.18 + /// \param _CapacityMap Type of capacity map.
10.19 + template <typename _Graph, typename _CapacityMap>
10.20 + struct EdmondsKarpDefaultTraits {
10.21 +
10.22 + /// \brief The graph type the algorithm runs on.
10.23 + typedef _Graph Graph;
10.24 +
10.25 + /// \brief The type of the map that stores the edge capacities.
10.26 + ///
10.27 + /// The type of the map that stores the edge capacities.
10.28 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
10.29 + typedef _CapacityMap CapacityMap;
10.30 +
10.31 + /// \brief The type of the length of the edges.
10.32 + typedef typename CapacityMap::Value Value;
10.33 +
10.34 + /// \brief The map type that stores the flow values.
10.35 + ///
10.36 + /// The map type that stores the flow values.
10.37 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
10.38 + typedef typename Graph::template EdgeMap<Value> FlowMap;
10.39 +
10.40 + /// \brief Instantiates a FlowMap.
10.41 + ///
10.42 + /// This function instantiates a \ref FlowMap.
10.43 + /// \param graph The graph, to which we would like to define the flow map.
10.44 + static FlowMap* createFlowMap(const Graph& graph) {
10.45 + return new FlowMap(graph);
10.46 + }
10.47 +
10.48 + /// \brief The tolerance used by the algorithm
10.49 + ///
10.50 + /// The tolerance used by the algorithm to handle inexact computation.
10.51 + typedef Tolerance<Value> Tolerance;
10.52 +
10.53 + };
10.54 +
10.55 /// \ingroup max_flow
10.56 + ///
10.57 /// \brief Edmonds-Karp algorithms class.
10.58 ///
10.59 /// This class provides an implementation of the \e Edmonds-Karp \e
10.60 /// algorithm producing a flow of maximum value in a directed
10.61 - /// graph. The Edmonds-Karp algorithm is slower than the Preflow algorithm
10.62 - /// but it has an advantage of the step-by-step execution control with
10.63 - /// feasible flow solutions. The \e source node, the \e target node, the \e
10.64 - /// capacity of the edges and the \e starting \e flow value of the
10.65 - /// edges should be passed to the algorithm through the
10.66 - /// constructor.
10.67 + /// graphs. The Edmonds-Karp algorithm is slower than the Preflow
10.68 + /// algorithm but it has an advantage of the step-by-step execution
10.69 + /// control with feasible flow solutions. The \e source node, the \e
10.70 + /// target node, the \e capacity of the edges and the \e starting \e
10.71 + /// flow value of the edges should be passed to the algorithm
10.72 + /// through the constructor.
10.73 ///
10.74 - /// The time complexity of the algorithm is \f$ O(n * e^2) \f$ in
10.75 + /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
10.76 /// worst case. Always try the preflow algorithm instead of this if
10.77 /// you just want to compute the optimal flow.
10.78 ///
10.79 /// \param _Graph The directed graph type the algorithm runs on.
10.80 - /// \param _Number The number type of the capacities and the flow values.
10.81 /// \param _CapacityMap The capacity map type.
10.82 - /// \param _FlowMap The flow map type.
10.83 - /// \param _Tolerance The tolerance class to handle computation problems.
10.84 + /// \param _Traits Traits class to set various data types used by
10.85 + /// the algorithm. The default traits class is \ref
10.86 + /// EdmondsKarpDefaultTraits. See \ref EdmondsKarpDefaultTraits for the
10.87 + /// documentation of a Edmonds-Karp traits class.
10.88 ///
10.89 /// \author Balazs Dezso
10.90 #ifdef DOXYGEN
10.91 - template <typename _Graph, typename _Number,
10.92 - typename _CapacityMap, typename _FlowMap, typename _Tolerance>
10.93 -#else
10.94 - template <typename _Graph, typename _Number,
10.95 - typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,
10.96 - typename _FlowMap = typename _Graph::template EdgeMap<_Number>,
10.97 - typename _Tolerance = Tolerance<_Number> >
10.98 + template <typename _Graph, typename _CapacityMap, typename _Traits>
10.99 +#else
10.100 + template <typename _Graph,
10.101 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
10.102 + typename _Traits = EdmondsKarpDefaultTraits<_Graph, _CapacityMap> >
10.103 #endif
10.104 class EdmondsKarp {
10.105 public:
10.106
10.107 + typedef _Traits Traits;
10.108 + typedef typename Traits::Graph Graph;
10.109 + typedef typename Traits::CapacityMap CapacityMap;
10.110 + typedef typename Traits::Value Value;
10.111 +
10.112 + typedef typename Traits::FlowMap FlowMap;
10.113 + typedef typename Traits::Tolerance Tolerance;
10.114 +
10.115 /// \brief \ref Exception for the case when the source equals the target.
10.116 ///
10.117 /// \ref Exception for the case when the source equals the target.
10.118 @@ -76,110 +124,239 @@
10.119 };
10.120
10.121
10.122 - /// \brief The graph type the algorithm runs on.
10.123 - typedef _Graph Graph;
10.124 - /// \brief The value type of the algorithms.
10.125 - typedef _Number Number;
10.126 - /// \brief The capacity map on the edges.
10.127 - typedef _CapacityMap CapacityMap;
10.128 - /// \brief The flow map on the edges.
10.129 - typedef _FlowMap FlowMap;
10.130 - /// \brief The tolerance used by the algorithm.
10.131 - typedef _Tolerance Tolerance;
10.132 -
10.133 - typedef ResGraphAdaptor<Graph, Number, CapacityMap,
10.134 - FlowMap, Tolerance> ResGraph;
10.135 -
10.136 private:
10.137
10.138 - typedef typename Graph::Node Node;
10.139 - typedef typename Graph::Edge Edge;
10.140 + GRAPH_TYPEDEFS(typename Graph);
10.141 + typedef typename Graph::template NodeMap<Edge> PredMap;
10.142
10.143 - typedef typename Graph::NodeIt NodeIt;
10.144 - typedef typename Graph::EdgeIt EdgeIt;
10.145 - typedef typename Graph::InEdgeIt InEdgeIt;
10.146 - typedef typename Graph::OutEdgeIt OutEdgeIt;
10.147 + const Graph& _graph;
10.148 + const CapacityMap* _capacity;
10.149 +
10.150 + Node _source, _target;
10.151 +
10.152 + FlowMap* _flow;
10.153 + bool _local_flow;
10.154 +
10.155 + PredMap* _pred;
10.156 + std::vector<Node> _queue;
10.157 +
10.158 + Tolerance _tolerance;
10.159 + Value _flow_value;
10.160 +
10.161 + void createStructures() {
10.162 + if (!_flow) {
10.163 + _flow = Traits::createFlowMap(_graph);
10.164 + _local_flow = true;
10.165 + }
10.166 + if (!_pred) {
10.167 + _pred = new PredMap(_graph);
10.168 + }
10.169 + _queue.resize(countNodes(_graph));
10.170 + }
10.171 +
10.172 + void destroyStructures() {
10.173 + if (_local_flow) {
10.174 + delete _flow;
10.175 + }
10.176 + if (_pred) {
10.177 + delete _pred;
10.178 + }
10.179 + }
10.180
10.181 public:
10.182
10.183 + ///\name Named template parameters
10.184 +
10.185 + ///@{
10.186 +
10.187 + template <typename _FlowMap>
10.188 + struct DefFlowMapTraits : public Traits {
10.189 + typedef _FlowMap FlowMap;
10.190 + static FlowMap *createFlowMap(const Graph&) {
10.191 + throw UninitializedParameter();
10.192 + }
10.193 + };
10.194 +
10.195 + /// \brief \ref named-templ-param "Named parameter" for setting
10.196 + /// FlowMap type
10.197 + ///
10.198 + /// \ref named-templ-param "Named parameter" for setting FlowMap
10.199 + /// type
10.200 + template <typename _FlowMap>
10.201 + struct DefFlowMap
10.202 + : public EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
10.203 + typedef EdmondsKarp<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> >
10.204 + Create;
10.205 + };
10.206 +
10.207 +
10.208 + /// @}
10.209 +
10.210 /// \brief The constructor of the class.
10.211 ///
10.212 /// The constructor of the class.
10.213 /// \param graph The directed graph the algorithm runs on.
10.214 + /// \param capacity The capacity of the edges.
10.215 /// \param source The source node.
10.216 /// \param target The target node.
10.217 - /// \param capacity The capacity of the edges.
10.218 - /// \param flow The flow of the edges.
10.219 - /// \param tolerance Tolerance class.
10.220 - EdmondsKarp(const Graph& graph, Node source, Node target,
10.221 - const CapacityMap& capacity, FlowMap& flow,
10.222 - const Tolerance& tolerance = Tolerance())
10.223 - : _graph(graph), _capacity(capacity), _flow(flow),
10.224 - _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance),
10.225 - _source(source), _target(target)
10.226 + EdmondsKarp(const Graph& graph, const CapacityMap& capacity,
10.227 + Node source, Node target)
10.228 + : _graph(graph), _capacity(&capacity), _source(source), _target(target),
10.229 + _flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
10.230 {
10.231 if (_source == _target) {
10.232 throw InvalidArgument();
10.233 }
10.234 }
10.235
10.236 + /// \brief Destrcutor.
10.237 + ///
10.238 + /// Destructor.
10.239 + ~EdmondsKarp() {
10.240 + destroyStructures();
10.241 + }
10.242 +
10.243 + /// \brief Sets the capacity map.
10.244 + ///
10.245 + /// Sets the capacity map.
10.246 + /// \return \c (*this)
10.247 + EdmondsKarp& capacityMap(const CapacityMap& map) {
10.248 + _capacity = ↦
10.249 + return *this;
10.250 + }
10.251 +
10.252 + /// \brief Sets the flow map.
10.253 + ///
10.254 + /// Sets the flow map.
10.255 + /// \return \c (*this)
10.256 + EdmondsKarp& flowMap(FlowMap& map) {
10.257 + if (_local_flow) {
10.258 + delete _flow;
10.259 + _local_flow = false;
10.260 + }
10.261 + _flow = ↦
10.262 + return *this;
10.263 + }
10.264 +
10.265 + /// \brief Returns the flow map.
10.266 + ///
10.267 + /// \return The flow map.
10.268 + const FlowMap& flowMap() {
10.269 + return *_flow;
10.270 + }
10.271 +
10.272 + /// \brief Sets the source node.
10.273 + ///
10.274 + /// Sets the source node.
10.275 + /// \return \c (*this)
10.276 + EdmondsKarp& source(const Node& node) {
10.277 + _source = node;
10.278 + return *this;
10.279 + }
10.280 +
10.281 + /// \brief Sets the target node.
10.282 + ///
10.283 + /// Sets the target node.
10.284 + /// \return \c (*this)
10.285 + EdmondsKarp& target(const Node& node) {
10.286 + _target = node;
10.287 + return *this;
10.288 + }
10.289 +
10.290 + /// \brief Sets the tolerance used by algorithm.
10.291 + ///
10.292 + /// Sets the tolerance used by algorithm.
10.293 + EdmondsKarp& tolerance(const Tolerance& tolerance) const {
10.294 + _tolerance = tolerance;
10.295 + return *this;
10.296 + }
10.297 +
10.298 + /// \brief Returns the tolerance used by algorithm.
10.299 + ///
10.300 + /// Returns the tolerance used by algorithm.
10.301 + const Tolerance& tolerance() const {
10.302 + return tolerance;
10.303 + }
10.304 +
10.305 + /// \name Execution control The simplest way to execute the
10.306 + /// algorithm is to use the \c run() member functions.
10.307 + /// \n
10.308 + /// If you need more control on initial solution or
10.309 + /// execution then you have to call one \ref init() function and then
10.310 + /// the start() or multiple times the \c augment() member function.
10.311 +
10.312 + ///@{
10.313 +
10.314 /// \brief Initializes the algorithm
10.315 ///
10.316 /// It sets the flow to empty flow.
10.317 void init() {
10.318 + createStructures();
10.319 for (EdgeIt it(_graph); it != INVALID; ++it) {
10.320 - _flow.set(it, 0);
10.321 + _flow->set(it, 0);
10.322 }
10.323 - _value = 0;
10.324 + _flow_value = 0;
10.325 }
10.326
10.327 /// \brief Initializes the algorithm
10.328 ///
10.329 - /// If the flow map initially flow this let the flow map
10.330 - /// unchanged but the flow value will be set by the flow
10.331 - /// on the outedges from the source.
10.332 - void flowInit() {
10.333 - _value = 0;
10.334 + /// Initializes the flow to the \c flowMap. The \c flowMap should
10.335 + /// contain a feasible flow, ie. in each node excluding the source
10.336 + /// and the target the incoming flow should be equal to the
10.337 + /// outgoing flow.
10.338 + template <typename FlowMap>
10.339 + void flowInit(const FlowMap& flowMap) {
10.340 + createStructures();
10.341 + for (EdgeIt e(_graph); e != INVALID; ++e) {
10.342 + _flow->set(e, flowMap[e]);
10.343 + }
10.344 + _flow_value = 0;
10.345 for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.346 - _value += _flow[jt];
10.347 + _flow_value += (*_flow)[jt];
10.348 }
10.349 for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.350 - _value -= _flow[jt];
10.351 + _flow_value -= (*_flow)[jt];
10.352 }
10.353 }
10.354
10.355 /// \brief Initializes the algorithm
10.356 ///
10.357 - /// If the flow map initially flow this let the flow map
10.358 - /// unchanged but the flow value will be set by the flow
10.359 - /// on the outedges from the source. It also checks that
10.360 - /// the flow map really contains a flow.
10.361 - /// \return %True when the flow map really a flow.
10.362 - bool checkedFlowInit() {
10.363 - _value = 0;
10.364 - for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.365 - _value += _flow[jt];
10.366 - }
10.367 - for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.368 - _value -= _flow[jt];
10.369 + /// Initializes the flow to the \c flowMap. The \c flowMap should
10.370 + /// contain a feasible flow, ie. in each node excluding the source
10.371 + /// and the target the incoming flow should be equal to the
10.372 + /// outgoing flow.
10.373 + /// \return %False when the given flowMap does not contain
10.374 + /// feasible flow.
10.375 + template <typename FlowMap>
10.376 + bool checkedFlowInit(const FlowMap& flowMap) {
10.377 + createStructures();
10.378 + for (EdgeIt e(_graph); e != INVALID; ++e) {
10.379 + _flow->set(e, flowMap[e]);
10.380 }
10.381 for (NodeIt it(_graph); it != INVALID; ++it) {
10.382 if (it == _source || it == _target) continue;
10.383 - Number outFlow = 0;
10.384 + Value outFlow = 0;
10.385 for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
10.386 - outFlow += _flow[jt];
10.387 + outFlow += (*_flow)[jt];
10.388 }
10.389 - Number inFlow = 0;
10.390 + Value inFlow = 0;
10.391 for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
10.392 - inFlow += _flow[jt];
10.393 + inFlow += (*_flow)[jt];
10.394 }
10.395 if (_tolerance.different(outFlow, inFlow)) {
10.396 return false;
10.397 }
10.398 }
10.399 for (EdgeIt it(_graph); it != INVALID; ++it) {
10.400 - if (_tolerance.less(_flow[it], 0)) return false;
10.401 - if (_tolerance.less(_capacity[it], _flow[it])) return false;
10.402 + if (_tolerance.less((*_flow)[it], 0)) return false;
10.403 + if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
10.404 + }
10.405 + _flow_value = 0;
10.406 + for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.407 + _flow_value += (*_flow)[jt];
10.408 + }
10.409 + for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
10.410 + _flow_value -= (*_flow)[jt];
10.411 }
10.412 return true;
10.413 }
10.414 @@ -195,32 +372,76 @@
10.415 /// \return %False when the augmenting is not success so the
10.416 /// current flow is a feasible and optimal solution.
10.417 bool augment() {
10.418 - typename Bfs<ResGraph>
10.419 - ::template DefDistMap<NullMap<Node, int> >
10.420 - ::Create bfs(_resgraph);
10.421 + for (NodeIt n(_graph); n != INVALID; ++n) {
10.422 + _pred->set(n, INVALID);
10.423 + }
10.424 +
10.425 + int first = 0, last = 1;
10.426 +
10.427 + _queue[0] = _source;
10.428 + _pred->set(_source, OutEdgeIt(_graph, _source));
10.429
10.430 - NullMap<Node, int> distMap;
10.431 - bfs.distMap(distMap);
10.432 -
10.433 - bfs.init();
10.434 - bfs.addSource(_source);
10.435 - bfs.start(_target);
10.436 + while (first != last && (*_pred)[_target] == INVALID) {
10.437 + Node n = _queue[first++];
10.438 +
10.439 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
10.440 + Value rem = (*_capacity)[e] - (*_flow)[e];
10.441 + Node t = _graph.target(e);
10.442 + if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
10.443 + _pred->set(t, e);
10.444 + _queue[last++] = t;
10.445 + }
10.446 + }
10.447 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
10.448 + Value rem = (*_flow)[e];
10.449 + Node t = _graph.source(e);
10.450 + if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
10.451 + _pred->set(t, e);
10.452 + _queue[last++] = t;
10.453 + }
10.454 + }
10.455 + }
10.456
10.457 - if (!bfs.reached(_target)) {
10.458 - return false;
10.459 + if ((*_pred)[_target] != INVALID) {
10.460 + Node n = _target;
10.461 + Edge e = (*_pred)[n];
10.462 +
10.463 + Value prem = (*_capacity)[e] - (*_flow)[e];
10.464 + n = _graph.source(e);
10.465 + while (n != _source) {
10.466 + e = (*_pred)[n];
10.467 + if (_graph.target(e) == n) {
10.468 + Value rem = (*_capacity)[e] - (*_flow)[e];
10.469 + if (rem < prem) prem = rem;
10.470 + n = _graph.source(e);
10.471 + } else {
10.472 + Value rem = (*_flow)[e];
10.473 + if (rem < prem) prem = rem;
10.474 + n = _graph.target(e);
10.475 + }
10.476 + }
10.477 +
10.478 + n = _target;
10.479 + e = (*_pred)[n];
10.480 +
10.481 + _flow->set(e, (*_flow)[e] + prem);
10.482 + n = _graph.source(e);
10.483 + while (n != _source) {
10.484 + e = (*_pred)[n];
10.485 + if (_graph.target(e) == n) {
10.486 + _flow->set(e, (*_flow)[e] + prem);
10.487 + n = _graph.source(e);
10.488 + } else {
10.489 + _flow->set(e, (*_flow)[e] - prem);
10.490 + n = _graph.target(e);
10.491 + }
10.492 + }
10.493 +
10.494 + _flow_value += prem;
10.495 + return true;
10.496 + } else {
10.497 + return false;
10.498 }
10.499 - Number min = _resgraph.rescap(bfs.predEdge(_target));
10.500 - for (Node it = bfs.predNode(_target); it != _source;
10.501 - it = bfs.predNode(it)) {
10.502 - if (min > _resgraph.rescap(bfs.predEdge(it))) {
10.503 - min = _resgraph.rescap(bfs.predEdge(it));
10.504 - }
10.505 - }
10.506 - for (Node it = _target; it != _source; it = bfs.predNode(it)) {
10.507 - _resgraph.augment(bfs.predEdge(it), min);
10.508 - }
10.509 - _value += min;
10.510 - return true;
10.511 }
10.512
10.513 /// \brief Executes the algorithm
10.514 @@ -230,13 +451,6 @@
10.515 while (augment()) {}
10.516 }
10.517
10.518 - /// \brief Gives back the current flow value.
10.519 - ///
10.520 - /// Gives back the current flow _value.
10.521 - Number flowValue() const {
10.522 - return _value;
10.523 - }
10.524 -
10.525 /// \brief runs the algorithm.
10.526 ///
10.527 /// It is just a shorthand for:
10.528 @@ -250,119 +464,59 @@
10.529 start();
10.530 }
10.531
10.532 + /// @}
10.533 +
10.534 + /// \name Query Functions
10.535 + /// The result of the %Dijkstra algorithm can be obtained using these
10.536 + /// functions.\n
10.537 + /// Before the use of these functions,
10.538 + /// either run() or start() must be called.
10.539 +
10.540 + ///@{
10.541 +
10.542 + /// \brief Returns the value of the maximum flow.
10.543 + ///
10.544 + /// Returns the value of the maximum flow by returning the excess
10.545 + /// of the target node \c t. This value equals to the value of
10.546 + /// the maximum flow already after the first phase.
10.547 + Value flowValue() const {
10.548 + return _flow_value;
10.549 + }
10.550 +
10.551 +
10.552 + /// \brief Returns the flow on the edge.
10.553 + ///
10.554 + /// Sets the \c flowMap to the flow on the edges. This method can
10.555 + /// be called after the second phase of algorithm.
10.556 + Value flow(const Edge& edge) const {
10.557 + return (*_flow)[edge];
10.558 + }
10.559 +
10.560 + /// \brief Returns true when the node is on the source side of minimum cut.
10.561 + ///
10.562 +
10.563 + /// Returns true when the node is on the source side of minimum
10.564 + /// cut. This method can be called both after running \ref
10.565 + /// startFirstPhase() and \ref startSecondPhase().
10.566 + bool minCut(const Node& node) const {
10.567 + return (*_pred)[node] != INVALID;
10.568 + }
10.569 +
10.570 /// \brief Returns a minimum value cut.
10.571 ///
10.572 /// Sets \c cut to the characteristic vector of a minimum value cut
10.573 /// It simply calls the minMinCut member.
10.574 /// \retval cut Write node bool map.
10.575 template <typename CutMap>
10.576 - void minCut(CutMap& cut) const {
10.577 - minMinCut(cut);
10.578 - }
10.579 + void minCutMap(CutMap& cutMap) const {
10.580 + for (NodeIt n(_graph); n != INVALID; ++n) {
10.581 + cutMap.set(n, (*_pred)[n] != INVALID);
10.582 + }
10.583 + cutMap.set(_source, true);
10.584 + }
10.585
10.586 - /// \brief Returns the inclusionwise minimum of the minimum value cuts.
10.587 - ///
10.588 - /// Sets \c cut to the characteristic vector of the minimum value cut
10.589 - /// which is inclusionwise minimum. It is computed by processing a
10.590 - /// bfs from the source node \c source in the residual graph.
10.591 - /// \retval cut Write node bool map.
10.592 - template <typename CutMap>
10.593 - void minMinCut(CutMap& cut) const {
10.594 + /// @}
10.595
10.596 - typename Bfs<ResGraph>
10.597 - ::template DefDistMap<NullMap<Node, int> >
10.598 - ::template DefProcessedMap<CutMap>
10.599 - ::Create bfs(_resgraph);
10.600 -
10.601 - NullMap<Node, int> distMap;
10.602 - bfs.distMap(distMap);
10.603 -
10.604 - bfs.processedMap(cut);
10.605 -
10.606 - bfs.run(_source);
10.607 - }
10.608 -
10.609 - /// \brief Returns the inclusionwise minimum of the minimum value cuts.
10.610 - ///
10.611 - /// Sets \c cut to the characteristic vector of the minimum value cut
10.612 - /// which is inclusionwise minimum. It is computed by processing a
10.613 - /// bfs from the source node \c source in the residual graph.
10.614 - /// \retval cut Write node bool map.
10.615 - template <typename CutMap>
10.616 - void maxMinCut(CutMap& cut) const {
10.617 -
10.618 - typedef RevGraphAdaptor<const ResGraph> RevGraph;
10.619 -
10.620 - RevGraph revgraph(_resgraph);
10.621 -
10.622 - typename Bfs<RevGraph>
10.623 - ::template DefDistMap<NullMap<Node, int> >
10.624 - ::template DefPredMap<NullMap<Node, Edge> >
10.625 - ::template DefProcessedMap<NotWriteMap<CutMap> >
10.626 - ::Create bfs(revgraph);
10.627 -
10.628 - NullMap<Node, int> distMap;
10.629 - bfs.distMap(distMap);
10.630 -
10.631 - NullMap<Node, Edge> predMap;
10.632 - bfs.predMap(predMap);
10.633 -
10.634 - NotWriteMap<CutMap> notcut(cut);
10.635 - bfs.processedMap(notcut);
10.636 -
10.637 - bfs.run(_target);
10.638 - }
10.639 -
10.640 - /// \brief Returns the source node.
10.641 - ///
10.642 - /// Returns the source node.
10.643 - ///
10.644 - Node source() const {
10.645 - return _source;
10.646 - }
10.647 -
10.648 - /// \brief Returns the target node.
10.649 - ///
10.650 - /// Returns the target node.
10.651 - ///
10.652 - Node target() const {
10.653 - return _target;
10.654 - }
10.655 -
10.656 - /// \brief Returns a reference to capacity map.
10.657 - ///
10.658 - /// Returns a reference to capacity map.
10.659 - ///
10.660 - const CapacityMap &capacityMap() const {
10.661 - return *_capacity;
10.662 - }
10.663 -
10.664 - /// \brief Returns a reference to flow map.
10.665 - ///
10.666 - /// Returns a reference to flow map.
10.667 - ///
10.668 - const FlowMap &flowMap() const {
10.669 - return *_flow;
10.670 - }
10.671 -
10.672 - /// \brief Returns the tolerance used by algorithm.
10.673 - ///
10.674 - /// Returns the tolerance used by algorithm.
10.675 - const Tolerance& tolerance() const {
10.676 - return tolerance;
10.677 - }
10.678 -
10.679 - private:
10.680 -
10.681 - const Graph& _graph;
10.682 - const CapacityMap& _capacity;
10.683 - FlowMap& _flow;
10.684 - Tolerance _tolerance;
10.685 -
10.686 - ResGraph _resgraph;
10.687 - Node _source, _target;
10.688 - Number _value;
10.689 -
10.690 };
10.691
10.692 }
11.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
11.2 +++ b/lemon/goldberg_tarjan.h Sat Nov 17 20:58:11 2007 +0000
11.3 @@ -0,0 +1,1020 @@
11.4 +/* -*- C++ -*-
11.5 + *
11.6 + * This file is a part of LEMON, a generic C++ optimization library
11.7 + *
11.8 + * Copyright (C) 2003-2007
11.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
11.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
11.11 + *
11.12 + * Permission to use, modify and distribute this software is granted
11.13 + * provided that this copyright notice appears in all copies. For
11.14 + * precise terms see the accompanying LICENSE file.
11.15 + *
11.16 + * This software is provided "AS IS" with no warranty of any kind,
11.17 + * express or implied, and with no claim as to its suitability for any
11.18 + * purpose.
11.19 + *
11.20 + */
11.21 +
11.22 +#ifndef LEMON_GOLDBERG_TARJAN_H
11.23 +#define LEMON_GOLDBERG_TARJAN_H
11.24 +
11.25 +#include <vector>
11.26 +#include <queue>
11.27 +
11.28 +#include <lemon/error.h>
11.29 +#include <lemon/bits/invalid.h>
11.30 +#include <lemon/tolerance.h>
11.31 +#include <lemon/maps.h>
11.32 +#include <lemon/graph_utils.h>
11.33 +#include <lemon/dynamic_tree.h>
11.34 +#include <limits>
11.35 +
11.36 +/// \file
11.37 +/// \ingroup max_flow
11.38 +/// \brief Implementation of the preflow algorithm.
11.39 +
11.40 +namespace lemon {
11.41 +
11.42 + /// \brief Default traits class of GoldbergTarjan class.
11.43 + ///
11.44 + /// Default traits class of GoldbergTarjan class.
11.45 + /// \param _Graph Graph type.
11.46 + /// \param _CapacityMap Type of capacity map.
11.47 + template <typename _Graph, typename _CapacityMap>
11.48 + struct GoldbergTarjanDefaultTraits {
11.49 +
11.50 + /// \brief The graph type the algorithm runs on.
11.51 + typedef _Graph Graph;
11.52 +
11.53 + /// \brief The type of the map that stores the edge capacities.
11.54 + ///
11.55 + /// The type of the map that stores the edge capacities.
11.56 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
11.57 + typedef _CapacityMap CapacityMap;
11.58 +
11.59 + /// \brief The type of the length of the edges.
11.60 + typedef typename CapacityMap::Value Value;
11.61 +
11.62 + /// \brief The map type that stores the flow values.
11.63 + ///
11.64 + /// The map type that stores the flow values.
11.65 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
11.66 + typedef typename Graph::template EdgeMap<Value> FlowMap;
11.67 +
11.68 + /// \brief Instantiates a FlowMap.
11.69 + ///
11.70 + /// This function instantiates a \ref FlowMap.
11.71 + /// \param graph The graph, to which we would like to define the flow map.
11.72 + static FlowMap* createFlowMap(const Graph& graph) {
11.73 + return new FlowMap(graph);
11.74 + }
11.75 +
11.76 + /// \brief The eleavator type used by GoldbergTarjan algorithm.
11.77 + ///
11.78 + /// The elevator type used by GoldbergTarjan algorithm.
11.79 + ///
11.80 + /// \sa Elevator
11.81 + /// \sa LinkedElevator
11.82 + typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
11.83 +
11.84 + /// \brief Instantiates an Elevator.
11.85 + ///
11.86 + /// This function instantiates a \ref Elevator.
11.87 + /// \param graph The graph, to which we would like to define the elevator.
11.88 + /// \param max_level The maximum level of the elevator.
11.89 + static Elevator* createElevator(const Graph& graph, int max_level) {
11.90 + return new Elevator(graph, max_level);
11.91 + }
11.92 +
11.93 + /// \brief The tolerance used by the algorithm
11.94 + ///
11.95 + /// The tolerance used by the algorithm to handle inexact computation.
11.96 + typedef Tolerance<Value> Tolerance;
11.97 +
11.98 + };
11.99 +
11.100 + /// \ingroup max_flow
11.101 + /// \brief Goldberg-Tarjan algorithms class.
11.102 + ///
11.103 + /// This class provides an implementation of the \e GoldbergTarjan
11.104 + /// \e algorithm producing a flow of maximum value in a directed
11.105 + /// graph. The GoldbergTarjan algorithm is a theoretical improvement
11.106 + /// of the Goldberg's \ref Preflow "preflow" algorithm by using the \ref
11.107 + /// DynamicTree "dynamic tree" data structure of Sleator and Tarjan.
11.108 + ///
11.109 + /// The original preflow algorithm with \e "highest label" or \e
11.110 + /// FIFO heuristic has \f$O(n^3)\f$ time complexity. The current
11.111 + /// algorithm improved this complexity to
11.112 + /// \f$O(nm\log(\frac{n^2}{m}))\f$. The algorithm builds limited
11.113 + /// size dynamic trees on the residual graph, which can be used to
11.114 + /// make multilevel push operations and gives a better bound on the
11.115 + /// number of non-saturating pushes.
11.116 + ///
11.117 + /// \param Graph The directed graph type the algorithm runs on.
11.118 + /// \param CapacityMap The capacity map type.
11.119 + /// \param _Traits Traits class to set various data types used by
11.120 + /// the algorithm. The default traits class is \ref
11.121 + /// GoldbergTarjanDefaultTraits. See \ref
11.122 + /// GoldbergTarjanDefaultTraits for the documentation of a
11.123 + /// Goldberg-Tarjan traits class.
11.124 + ///
11.125 + /// \author Tamas Hamori and Balazs Dezso
11.126 +#ifdef DOXYGEN
11.127 + template <typename _Graph, typename _CapacityMap, typename _Traits>
11.128 +#else
11.129 + template <typename _Graph,
11.130 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
11.131 + typename _Traits =
11.132 + GoldbergTarjanDefaultTraits<_Graph, _CapacityMap> >
11.133 +#endif
11.134 + class GoldbergTarjan {
11.135 + public:
11.136 +
11.137 + typedef _Traits Traits;
11.138 + typedef typename Traits::Graph Graph;
11.139 + typedef typename Traits::CapacityMap CapacityMap;
11.140 + typedef typename Traits::Value Value;
11.141 +
11.142 + typedef typename Traits::FlowMap FlowMap;
11.143 + typedef typename Traits::Elevator Elevator;
11.144 + typedef typename Traits::Tolerance Tolerance;
11.145 +
11.146 + protected:
11.147 +
11.148 + GRAPH_TYPEDEFS(typename Graph);
11.149 +
11.150 + typedef typename Graph::template NodeMap<Node> NodeNodeMap;
11.151 + typedef typename Graph::template NodeMap<int> IntNodeMap;
11.152 +
11.153 + typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
11.154 + typedef typename Graph::template EdgeMap<Edge> EdgeEdgeMap;
11.155 +
11.156 + typedef typename std::vector<Node> VecNode;
11.157 +
11.158 + typedef DynamicTree<Value,IntNodeMap,Tolerance> DynTree;
11.159 +
11.160 + const Graph& _graph;
11.161 + const CapacityMap* _capacity;
11.162 + int _node_num; //the number of nodes of G
11.163 +
11.164 + Node _source;
11.165 + Node _target;
11.166 +
11.167 + FlowMap* _flow;
11.168 + bool _local_flow;
11.169 +
11.170 + Elevator* _level;
11.171 + bool _local_level;
11.172 +
11.173 + typedef typename Graph::template NodeMap<Value> ExcessMap;
11.174 + ExcessMap* _excess;
11.175 +
11.176 + Tolerance _tolerance;
11.177 +
11.178 + bool _phase;
11.179 +
11.180 + // constant for treesize
11.181 + static const double _tree_bound = 2;
11.182 + int _max_tree_size;
11.183 +
11.184 + //tags for the dynamic tree
11.185 + DynTree *_dt;
11.186 + //datastructure of dyanamic tree
11.187 + IntNodeMap *_dt_index;
11.188 + //datastrucure for solution of the communication between the two class
11.189 + EdgeNodeMap *_dt_edges;
11.190 + //nodeMap for storing the outgoing edge from the node in the tree
11.191 +
11.192 + //max of the type Value
11.193 + const Value _max_value;
11.194 +
11.195 + public:
11.196 +
11.197 + typedef GoldbergTarjan Create;
11.198 +
11.199 + ///\name Named template parameters
11.200 +
11.201 + ///@{
11.202 +
11.203 + template <typename _FlowMap>
11.204 + struct DefFlowMapTraits : public Traits {
11.205 + typedef _FlowMap FlowMap;
11.206 + static FlowMap *createFlowMap(const Graph&) {
11.207 + throw UninitializedParameter();
11.208 + }
11.209 + };
11.210 +
11.211 + /// \brief \ref named-templ-param "Named parameter" for setting
11.212 + /// FlowMap type
11.213 + ///
11.214 + /// \ref named-templ-param "Named parameter" for setting FlowMap
11.215 + /// type
11.216 + template <typename _FlowMap>
11.217 + struct DefFlowMap
11.218 + : public GoldbergTarjan<Graph, CapacityMap,
11.219 + DefFlowMapTraits<_FlowMap> > {
11.220 + typedef GoldbergTarjan<Graph, CapacityMap,
11.221 + DefFlowMapTraits<_FlowMap> > Create;
11.222 + };
11.223 +
11.224 + template <typename _Elevator>
11.225 + struct DefElevatorTraits : public Traits {
11.226 + typedef _Elevator Elevator;
11.227 + static Elevator *createElevator(const Graph&, int) {
11.228 + throw UninitializedParameter();
11.229 + }
11.230 + };
11.231 +
11.232 + /// \brief \ref named-templ-param "Named parameter" for setting
11.233 + /// Elevator type
11.234 + ///
11.235 + /// \ref named-templ-param "Named parameter" for setting Elevator
11.236 + /// type
11.237 + template <typename _Elevator>
11.238 + struct DefElevator
11.239 + : public GoldbergTarjan<Graph, CapacityMap,
11.240 + DefElevatorTraits<_Elevator> > {
11.241 + typedef GoldbergTarjan<Graph, CapacityMap,
11.242 + DefElevatorTraits<_Elevator> > Create;
11.243 + };
11.244 +
11.245 + template <typename _Elevator>
11.246 + struct DefStandardElevatorTraits : public Traits {
11.247 + typedef _Elevator Elevator;
11.248 + static Elevator *createElevator(const Graph& graph, int max_level) {
11.249 + return new Elevator(graph, max_level);
11.250 + }
11.251 + };
11.252 +
11.253 + /// \brief \ref named-templ-param "Named parameter" for setting
11.254 + /// Elevator type
11.255 + ///
11.256 + /// \ref named-templ-param "Named parameter" for setting Elevator
11.257 + /// type. The Elevator should be standard constructor interface, ie.
11.258 + /// the graph and the maximum level should be passed to it.
11.259 + template <typename _Elevator>
11.260 + struct DefStandardElevator
11.261 + : public GoldbergTarjan<Graph, CapacityMap,
11.262 + DefStandardElevatorTraits<_Elevator> > {
11.263 + typedef GoldbergTarjan<Graph, CapacityMap,
11.264 + DefStandardElevatorTraits<_Elevator> > Create;
11.265 + };
11.266 +
11.267 +
11.268 + ///\ref Exception for the case when s=t.
11.269 +
11.270 + ///\ref Exception for the case when the source equals the target.
11.271 + class InvalidArgument : public lemon::LogicError {
11.272 + public:
11.273 + virtual const char* what() const throw() {
11.274 + return "lemon::GoldbergTarjan::InvalidArgument";
11.275 + }
11.276 + };
11.277 +
11.278 + public:
11.279 +
11.280 + /// \brief The constructor of the class.
11.281 + ///
11.282 + /// The constructor of the class.
11.283 + /// \param graph The directed graph the algorithm runs on.
11.284 + /// \param capacity The capacity of the edges.
11.285 + /// \param source The source node.
11.286 + /// \param target The target node.
11.287 + /// Except the graph, all of these parameters can be reset by
11.288 + /// calling \ref source, \ref target, \ref capacityMap and \ref
11.289 + /// flowMap, resp.
11.290 + GoldbergTarjan(const Graph& graph, const CapacityMap& capacity,
11.291 + Node source, Node target)
11.292 + : _graph(graph), _capacity(&capacity),
11.293 + _node_num(), _source(source), _target(target),
11.294 + _flow(0), _local_flow(false),
11.295 + _level(0), _local_level(false),
11.296 + _excess(0), _tolerance(),
11.297 + _phase(), _max_tree_size(),
11.298 + _dt(0), _dt_index(0), _dt_edges(0),
11.299 + _max_value(std::numeric_limits<Value>::max()) {
11.300 + if (_source == _target) throw InvalidArgument();
11.301 + }
11.302 +
11.303 + /// \brief Destrcutor.
11.304 + ///
11.305 + /// Destructor.
11.306 + ~GoldbergTarjan() {
11.307 + destroyStructures();
11.308 + }
11.309 +
11.310 + /// \brief Sets the capacity map.
11.311 + ///
11.312 + /// Sets the capacity map.
11.313 + /// \return \c (*this)
11.314 + GoldbergTarjan& capacityMap(const CapacityMap& map) {
11.315 + _capacity = ↦
11.316 + return *this;
11.317 + }
11.318 +
11.319 + /// \brief Sets the flow map.
11.320 + ///
11.321 + /// Sets the flow map.
11.322 + /// \return \c (*this)
11.323 + GoldbergTarjan& flowMap(FlowMap& map) {
11.324 + if (_local_flow) {
11.325 + delete _flow;
11.326 + _local_flow = false;
11.327 + }
11.328 + _flow = ↦
11.329 + return *this;
11.330 + }
11.331 +
11.332 + /// \brief Returns the flow map.
11.333 + ///
11.334 + /// \return The flow map.
11.335 + const FlowMap& flowMap() {
11.336 + return *_flow;
11.337 + }
11.338 +
11.339 + /// \brief Sets the elevator.
11.340 + ///
11.341 + /// Sets the elevator.
11.342 + /// \return \c (*this)
11.343 + GoldbergTarjan& elevator(Elevator& elevator) {
11.344 + if (_local_level) {
11.345 + delete _level;
11.346 + _local_level = false;
11.347 + }
11.348 + _level = &elevator;
11.349 + return *this;
11.350 + }
11.351 +
11.352 + /// \brief Returns the elevator.
11.353 + ///
11.354 + /// \return The elevator.
11.355 + const Elevator& elevator() {
11.356 + return *_level;
11.357 + }
11.358 +
11.359 + /// \brief Sets the source node.
11.360 + ///
11.361 + /// Sets the source node.
11.362 + /// \return \c (*this)
11.363 + GoldbergTarjan& source(const Node& node) {
11.364 + _source = node;
11.365 + return *this;
11.366 + }
11.367 +
11.368 + /// \brief Sets the target node.
11.369 + ///
11.370 + /// Sets the target node.
11.371 + /// \return \c (*this)
11.372 + GoldbergTarjan& target(const Node& node) {
11.373 + _target = node;
11.374 + return *this;
11.375 + }
11.376 +
11.377 + /// \brief Sets the tolerance used by algorithm.
11.378 + ///
11.379 + /// Sets the tolerance used by algorithm.
11.380 + GoldbergTarjan& tolerance(const Tolerance& tolerance) const {
11.381 + _tolerance = tolerance;
11.382 + if (_dt) {
11.383 + _dt->tolerance(_tolerance);
11.384 + }
11.385 + return *this;
11.386 + }
11.387 +
11.388 + /// \brief Returns the tolerance used by algorithm.
11.389 + ///
11.390 + /// Returns the tolerance used by algorithm.
11.391 + const Tolerance& tolerance() const {
11.392 + return tolerance;
11.393 + }
11.394 +
11.395 +
11.396 + private:
11.397 +
11.398 + void createStructures() {
11.399 + _node_num = countNodes(_graph);
11.400 +
11.401 + _max_tree_size = (double(_node_num) * double(_node_num)) /
11.402 + double(countEdges(_graph));
11.403 +
11.404 + if (!_flow) {
11.405 + _flow = Traits::createFlowMap(_graph);
11.406 + _local_flow = true;
11.407 + }
11.408 + if (!_level) {
11.409 + _level = Traits::createElevator(_graph, _node_num);
11.410 + _local_level = true;
11.411 + }
11.412 + if (!_excess) {
11.413 + _excess = new ExcessMap(_graph);
11.414 + }
11.415 + if (!_dt_index && !_dt) {
11.416 + _dt_index = new IntNodeMap(_graph);
11.417 + _dt = new DynTree(*_dt_index, _tolerance);
11.418 + }
11.419 + if (!_dt_edges) {
11.420 + _dt_edges = new EdgeNodeMap(_graph);
11.421 + }
11.422 + }
11.423 +
11.424 + void destroyStructures() {
11.425 + if (_local_flow) {
11.426 + delete _flow;
11.427 + }
11.428 + if (_local_level) {
11.429 + delete _level;
11.430 + }
11.431 + if (_excess) {
11.432 + delete _excess;
11.433 + }
11.434 + if (_dt) {
11.435 + delete _dt;
11.436 + }
11.437 + if (_dt_index) {
11.438 + delete _dt_index;
11.439 + }
11.440 + if (_dt_edges) {
11.441 + delete _dt_edges;
11.442 + }
11.443 + }
11.444 +
11.445 + bool sendOut(Node n, Value& excess) {
11.446 +
11.447 + Node w = _dt->findRoot(n);
11.448 +
11.449 + while (w != n) {
11.450 +
11.451 + Value rem;
11.452 + Node u = _dt->findCost(n, rem);
11.453 +
11.454 + if (_tolerance.positive(rem) && !_level->active(w) && w != _target) {
11.455 + _level->activate(w);
11.456 + }
11.457 +
11.458 + if (!_tolerance.less(rem, excess)) {
11.459 + _dt->addCost(n, - excess);
11.460 + _excess->set(w, (*_excess)[w] + excess);
11.461 + excess = 0;
11.462 + return true;
11.463 + }
11.464 +
11.465 +
11.466 + _dt->addCost(n, - rem);
11.467 +
11.468 + excess -= rem;
11.469 + _excess->set(w, (*_excess)[w] + rem);
11.470 +
11.471 + _dt->cut(u);
11.472 + _dt->addCost(u, _max_value);
11.473 +
11.474 + Edge e = (*_dt_edges)[u];
11.475 + _dt_edges->set(u, INVALID);
11.476 +
11.477 + if (u == _graph.source(e)) {
11.478 + _flow->set(e, (*_capacity)[e]);
11.479 + } else {
11.480 + _flow->set(e, 0);
11.481 + }
11.482 +
11.483 + w = _dt->findRoot(n);
11.484 + }
11.485 + return false;
11.486 + }
11.487 +
11.488 + bool sendIn(Node n, Value& excess) {
11.489 +
11.490 + Node w = _dt->findRoot(n);
11.491 +
11.492 + while (w != n) {
11.493 +
11.494 + Value rem;
11.495 + Node u = _dt->findCost(n, rem);
11.496 +
11.497 + if (_tolerance.positive(rem) && !_level->active(w) && w != _source) {
11.498 + _level->activate(w);
11.499 + }
11.500 +
11.501 + if (!_tolerance.less(rem, excess)) {
11.502 + _dt->addCost(n, - excess);
11.503 + _excess->set(w, (*_excess)[w] + excess);
11.504 + excess = 0;
11.505 + return true;
11.506 + }
11.507 +
11.508 +
11.509 + _dt->addCost(n, - rem);
11.510 +
11.511 + excess -= rem;
11.512 + _excess->set(w, (*_excess)[w] + rem);
11.513 +
11.514 + _dt->cut(u);
11.515 + _dt->addCost(u, _max_value);
11.516 +
11.517 + Edge e = (*_dt_edges)[u];
11.518 + _dt_edges->set(u, INVALID);
11.519 +
11.520 + if (u == _graph.source(e)) {
11.521 + _flow->set(e, (*_capacity)[e]);
11.522 + } else {
11.523 + _flow->set(e, 0);
11.524 + }
11.525 +
11.526 + w = _dt->findRoot(n);
11.527 + }
11.528 + return false;
11.529 + }
11.530 +
11.531 + void cutChildren(Node n) {
11.532 +
11.533 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
11.534 +
11.535 + Node v = _graph.target(e);
11.536 +
11.537 + if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
11.538 + _dt->cut(v);
11.539 + _dt_edges->set(v, INVALID);
11.540 + Value rem;
11.541 + _dt->findCost(v, rem);
11.542 + _dt->addCost(v, - rem);
11.543 + _dt->addCost(v, _max_value);
11.544 + _flow->set(e, rem);
11.545 +
11.546 + }
11.547 + }
11.548 +
11.549 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
11.550 +
11.551 + Node v = _graph.source(e);
11.552 +
11.553 + if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
11.554 + _dt->cut(v);
11.555 + _dt_edges->set(v, INVALID);
11.556 + Value rem;
11.557 + _dt->findCost(v, rem);
11.558 + _dt->addCost(v, - rem);
11.559 + _dt->addCost(v, _max_value);
11.560 + _flow->set(e, (*_capacity)[e] - rem);
11.561 +
11.562 + }
11.563 + }
11.564 + }
11.565 +
11.566 + void extractTrees() {
11.567 + for (NodeIt n(_graph); n != INVALID; ++n) {
11.568 +
11.569 + Node w = _dt->findRoot(n);
11.570 +
11.571 + while (w != n) {
11.572 +
11.573 + Value rem;
11.574 + Node u = _dt->findCost(n, rem);
11.575 +
11.576 + _dt->cut(u);
11.577 + _dt->addCost(u, - rem);
11.578 + _dt->addCost(u, _max_value);
11.579 +
11.580 + Edge e = (*_dt_edges)[u];
11.581 + _dt_edges->set(u, INVALID);
11.582 +
11.583 + if (u == _graph.source(e)) {
11.584 + _flow->set(e, (*_capacity)[e] - rem);
11.585 + } else {
11.586 + _flow->set(e, rem);
11.587 + }
11.588 +
11.589 + w = _dt->findRoot(n);
11.590 + }
11.591 + }
11.592 + }
11.593 +
11.594 + public:
11.595 +
11.596 + /// \name Execution control The simplest way to execute the
11.597 + /// algorithm is to use one of the member functions called \c
11.598 + /// run().
11.599 + /// \n
11.600 + /// If you need more control on initial solution or
11.601 + /// execution then you have to call one \ref init() function and then
11.602 + /// the startFirstPhase() and if you need the startSecondPhase().
11.603 +
11.604 + ///@{
11.605 +
11.606 + /// \brief Initializes the internal data structures.
11.607 + ///
11.608 + /// Initializes the internal data structures.
11.609 + ///
11.610 + void init() {
11.611 + createStructures();
11.612 +
11.613 + for (NodeIt n(_graph); n != INVALID; ++n) {
11.614 + _excess->set(n, 0);
11.615 + }
11.616 +
11.617 + for (EdgeIt e(_graph); e != INVALID; ++e) {
11.618 + _flow->set(e, 0);
11.619 + }
11.620 +
11.621 + _dt->clear();
11.622 + for (NodeIt v(_graph); v != INVALID; ++v) {
11.623 + (*_dt_edges)[v] = INVALID;
11.624 + _dt->makeTree(v);
11.625 + _dt->addCost(v, _max_value);
11.626 + }
11.627 +
11.628 + typename Graph::template NodeMap<bool> reached(_graph, false);
11.629 +
11.630 + _level->initStart();
11.631 + _level->initAddItem(_target);
11.632 +
11.633 + std::vector<Node> queue;
11.634 + reached.set(_source, true);
11.635 +
11.636 + queue.push_back(_target);
11.637 + reached.set(_target, true);
11.638 + while (!queue.empty()) {
11.639 + _level->initNewLevel();
11.640 + std::vector<Node> nqueue;
11.641 + for (int i = 0; i < int(queue.size()); ++i) {
11.642 + Node n = queue[i];
11.643 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
11.644 + Node u = _graph.source(e);
11.645 + if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
11.646 + reached.set(u, true);
11.647 + _level->initAddItem(u);
11.648 + nqueue.push_back(u);
11.649 + }
11.650 + }
11.651 + }
11.652 + queue.swap(nqueue);
11.653 + }
11.654 + _level->initFinish();
11.655 +
11.656 + for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
11.657 + if (_tolerance.positive((*_capacity)[e])) {
11.658 + Node u = _graph.target(e);
11.659 + if ((*_level)[u] == _level->maxLevel()) continue;
11.660 + _flow->set(e, (*_capacity)[e]);
11.661 + _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
11.662 + if (u != _target && !_level->active(u)) {
11.663 + _level->activate(u);
11.664 + }
11.665 + }
11.666 + }
11.667 + }
11.668 +
11.669 + /// \brief Starts the first phase of the preflow algorithm.
11.670 + ///
11.671 + /// The preflow algorithm consists of two phases, this method runs
11.672 + /// the first phase. After the first phase the maximum flow value
11.673 + /// and a minimum value cut can already be computed, although a
11.674 + /// maximum flow is not yet obtained. So after calling this method
11.675 + /// \ref flowValue() returns the value of a maximum flow and \ref
11.676 + /// minCut() returns a minimum cut.
11.677 + /// \pre One of the \ref init() functions should be called.
11.678 + void startFirstPhase() {
11.679 + _phase = true;
11.680 + Node n;
11.681 +
11.682 + while ((n = _level->highestActive()) != INVALID) {
11.683 + Value excess = (*_excess)[n];
11.684 + int level = _level->highestActiveLevel();
11.685 + int new_level = _level->maxLevel();
11.686 +
11.687 + if (_dt->findRoot(n) != n) {
11.688 + if (sendOut(n, excess)) goto no_more_push;
11.689 + }
11.690 +
11.691 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
11.692 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.693 + Node v = _graph.target(e);
11.694 +
11.695 + if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
11.696 +
11.697 + if ((*_level)[v] < level) {
11.698 +
11.699 + if (_dt->findSize(n) + _dt->findSize(v) <
11.700 + _tree_bound * _max_tree_size) {
11.701 + _dt->addCost(n, -_max_value);
11.702 + _dt->addCost(n, rem);
11.703 + _dt->link(n, v);
11.704 + _dt_edges->set(n, e);
11.705 + if (sendOut(n, excess)) goto no_more_push;
11.706 + } else {
11.707 + if (!_level->active(v) && v != _target) {
11.708 + _level->activate(v);
11.709 + }
11.710 + if (!_tolerance.less(rem, excess)) {
11.711 + _flow->set(e, (*_flow)[e] + excess);
11.712 + _excess->set(v, (*_excess)[v] + excess);
11.713 + excess = 0;
11.714 + goto no_more_push;
11.715 + } else {
11.716 + excess -= rem;
11.717 + _excess->set(v, (*_excess)[v] + rem);
11.718 + _flow->set(e, (*_capacity)[e]);
11.719 + }
11.720 + }
11.721 + } else if (new_level > (*_level)[v]) {
11.722 + new_level = (*_level)[v];
11.723 + }
11.724 + }
11.725 +
11.726 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
11.727 + Value rem = (*_flow)[e];
11.728 + Node v = _graph.source(e);
11.729 + if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
11.730 +
11.731 + if ((*_level)[v] < level) {
11.732 +
11.733 + if (_dt->findSize(n) + _dt->findSize(v) <
11.734 + _tree_bound * _max_tree_size) {
11.735 + _dt->addCost(n, - _max_value);
11.736 + _dt->addCost(n, rem);
11.737 + _dt->link(n, v);
11.738 + _dt_edges->set(n, e);
11.739 + if (sendOut(n, excess)) goto no_more_push;
11.740 + } else {
11.741 + if (!_level->active(v) && v != _target) {
11.742 + _level->activate(v);
11.743 + }
11.744 + if (!_tolerance.less(rem, excess)) {
11.745 + _flow->set(e, (*_flow)[e] - excess);
11.746 + _excess->set(v, (*_excess)[v] + excess);
11.747 + excess = 0;
11.748 + goto no_more_push;
11.749 + } else {
11.750 + excess -= rem;
11.751 + _excess->set(v, (*_excess)[v] + rem);
11.752 + _flow->set(e, 0);
11.753 + }
11.754 + }
11.755 + } else if (new_level > (*_level)[v]) {
11.756 + new_level = (*_level)[v];
11.757 + }
11.758 + }
11.759 +
11.760 + no_more_push:
11.761 +
11.762 + _excess->set(n, excess);
11.763 +
11.764 + if (excess != 0) {
11.765 + cutChildren(n);
11.766 + if (new_level + 1 < _level->maxLevel()) {
11.767 + _level->liftHighestActive(new_level + 1);
11.768 + } else {
11.769 + _level->liftHighestActiveToTop();
11.770 + }
11.771 + if (_level->emptyLevel(level)) {
11.772 + _level->liftToTop(level);
11.773 + }
11.774 + } else {
11.775 + _level->deactivate(n);
11.776 + }
11.777 + }
11.778 + extractTrees();
11.779 + }
11.780 +
11.781 + /// \brief Starts the second phase of the preflow algorithm.
11.782 + ///
11.783 + /// The preflow algorithm consists of two phases, this method runs
11.784 + /// the second phase. After calling \ref init() and \ref
11.785 + /// startFirstPhase() and then \ref startSecondPhase(), \ref
11.786 + /// flowMap() return a maximum flow, \ref flowValue() returns the
11.787 + /// value of a maximum flow, \ref minCut() returns a minimum cut
11.788 + /// \pre The \ref init() and startFirstPhase() functions should be
11.789 + /// called before.
11.790 + void startSecondPhase() {
11.791 + _phase = false;
11.792 +
11.793 + typename Graph::template NodeMap<bool> reached(_graph, false);
11.794 + for (NodeIt n(_graph); n != INVALID; ++n) {
11.795 + reached.set(n, (*_level)[n] < _level->maxLevel());
11.796 + }
11.797 +
11.798 + _level->initStart();
11.799 + _level->initAddItem(_source);
11.800 +
11.801 + std::vector<Node> queue;
11.802 + queue.push_back(_source);
11.803 + reached.set(_source, true);
11.804 +
11.805 + while (!queue.empty()) {
11.806 + _level->initNewLevel();
11.807 + std::vector<Node> nqueue;
11.808 + for (int i = 0; i < int(queue.size()); ++i) {
11.809 + Node n = queue[i];
11.810 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
11.811 + Node v = _graph.target(e);
11.812 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
11.813 + reached.set(v, true);
11.814 + _level->initAddItem(v);
11.815 + nqueue.push_back(v);
11.816 + }
11.817 + }
11.818 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
11.819 + Node u = _graph.source(e);
11.820 + if (!reached[u] &&
11.821 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
11.822 + reached.set(u, true);
11.823 + _level->initAddItem(u);
11.824 + nqueue.push_back(u);
11.825 + }
11.826 + }
11.827 + }
11.828 + queue.swap(nqueue);
11.829 + }
11.830 + _level->initFinish();
11.831 +
11.832 + for (NodeIt n(_graph); n != INVALID; ++n) {
11.833 + if ((*_excess)[n] > 0 && _target != n) {
11.834 + _level->activate(n);
11.835 + }
11.836 + }
11.837 +
11.838 + Node n;
11.839 +
11.840 + while ((n = _level->highestActive()) != INVALID) {
11.841 + Value excess = (*_excess)[n];
11.842 + int level = _level->highestActiveLevel();
11.843 + int new_level = _level->maxLevel();
11.844 +
11.845 + if (_dt->findRoot(n) != n) {
11.846 + if (sendIn(n, excess)) goto no_more_push;
11.847 + }
11.848 +
11.849 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
11.850 + Value rem = (*_capacity)[e] - (*_flow)[e];
11.851 + Node v = _graph.target(e);
11.852 +
11.853 + if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
11.854 +
11.855 + if ((*_level)[v] < level) {
11.856 +
11.857 + if (_dt->findSize(n) + _dt->findSize(v) <
11.858 + _tree_bound * _max_tree_size) {
11.859 + _dt->addCost(n, -_max_value);
11.860 + _dt->addCost(n, rem);
11.861 + _dt->link(n, v);
11.862 + _dt_edges->set(n, e);
11.863 + if (sendIn(n, excess)) goto no_more_push;
11.864 + } else {
11.865 + if (!_level->active(v) && v != _source) {
11.866 + _level->activate(v);
11.867 + }
11.868 + if (!_tolerance.less(rem, excess)) {
11.869 + _flow->set(e, (*_flow)[e] + excess);
11.870 + _excess->set(v, (*_excess)[v] + excess);
11.871 + excess = 0;
11.872 + goto no_more_push;
11.873 + } else {
11.874 + excess -= rem;
11.875 + _excess->set(v, (*_excess)[v] + rem);
11.876 + _flow->set(e, (*_capacity)[e]);
11.877 + }
11.878 + }
11.879 + } else if (new_level > (*_level)[v]) {
11.880 + new_level = (*_level)[v];
11.881 + }
11.882 + }
11.883 +
11.884 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
11.885 + Value rem = (*_flow)[e];
11.886 + Node v = _graph.source(e);
11.887 + if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
11.888 +
11.889 + if ((*_level)[v] < level) {
11.890 +
11.891 + if (_dt->findSize(n) + _dt->findSize(v) <
11.892 + _tree_bound * _max_tree_size) {
11.893 + _dt->addCost(n, - _max_value);
11.894 + _dt->addCost(n, rem);
11.895 + _dt->link(n, v);
11.896 + _dt_edges->set(n, e);
11.897 + if (sendIn(n, excess)) goto no_more_push;
11.898 + } else {
11.899 + if (!_level->active(v) && v != _source) {
11.900 + _level->activate(v);
11.901 + }
11.902 + if (!_tolerance.less(rem, excess)) {
11.903 + _flow->set(e, (*_flow)[e] - excess);
11.904 + _excess->set(v, (*_excess)[v] + excess);
11.905 + excess = 0;
11.906 + goto no_more_push;
11.907 + } else {
11.908 + excess -= rem;
11.909 + _excess->set(v, (*_excess)[v] + rem);
11.910 + _flow->set(e, 0);
11.911 + }
11.912 + }
11.913 + } else if (new_level > (*_level)[v]) {
11.914 + new_level = (*_level)[v];
11.915 + }
11.916 + }
11.917 +
11.918 + no_more_push:
11.919 +
11.920 + _excess->set(n, excess);
11.921 +
11.922 + if (excess != 0) {
11.923 + cutChildren(n);
11.924 + if (new_level + 1 < _level->maxLevel()) {
11.925 + _level->liftHighestActive(new_level + 1);
11.926 + } else {
11.927 + _level->liftHighestActiveToTop();
11.928 + }
11.929 + if (_level->emptyLevel(level)) {
11.930 + _level->liftToTop(level);
11.931 + }
11.932 + } else {
11.933 + _level->deactivate(n);
11.934 + }
11.935 + }
11.936 + extractTrees();
11.937 + }
11.938 +
11.939 + /// \brief Runs the Goldberg-Tarjan algorithm.
11.940 + ///
11.941 + /// Runs the Goldberg-Tarjan algorithm.
11.942 + /// \note pf.run() is just a shortcut of the following code.
11.943 + /// \code
11.944 + /// pf.init();
11.945 + /// pf.startFirstPhase();
11.946 + /// pf.startSecondPhase();
11.947 + /// \endcode
11.948 + void run() {
11.949 + init();
11.950 + startFirstPhase();
11.951 + startSecondPhase();
11.952 + }
11.953 +
11.954 + /// \brief Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
11.955 + ///
11.956 + /// Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
11.957 + /// \note pf.runMinCut() is just a shortcut of the following code.
11.958 + /// \code
11.959 + /// pf.init();
11.960 + /// pf.startFirstPhase();
11.961 + /// \endcode
11.962 + void runMinCut() {
11.963 + init();
11.964 + startFirstPhase();
11.965 + }
11.966 +
11.967 + /// @}
11.968 +
11.969 + /// \name Query Functions
11.970 + /// The result of the %Dijkstra algorithm can be obtained using these
11.971 + /// functions.\n
11.972 + /// Before the use of these functions,
11.973 + /// either run() or start() must be called.
11.974 +
11.975 + ///@{
11.976 +
11.977 + /// \brief Returns the value of the maximum flow.
11.978 + ///
11.979 + /// Returns the value of the maximum flow by returning the excess
11.980 + /// of the target node \c t. This value equals to the value of
11.981 + /// the maximum flow already after the first phase.
11.982 + Value flowValue() const {
11.983 + return (*_excess)[_target];
11.984 + }
11.985 +
11.986 + /// \brief Returns true when the node is on the source side of minimum cut.
11.987 + ///
11.988 + /// Returns true when the node is on the source side of minimum
11.989 + /// cut. This method can be called both after running \ref
11.990 + /// startFirstPhase() and \ref startSecondPhase().
11.991 + bool minCut(const Node& node) const {
11.992 + return ((*_level)[node] == _level->maxLevel()) == _phase;
11.993 + }
11.994 +
11.995 + /// \brief Returns a minimum value cut.
11.996 + ///
11.997 + /// Sets the \c cutMap to the characteristic vector of a minimum value
11.998 + /// cut. This method can be called both after running \ref
11.999 + /// startFirstPhase() and \ref startSecondPhase(). The result after second
11.1000 + /// phase could be changed slightly if inexact computation is used.
11.1001 + /// \pre The \c cutMap should be a bool-valued node-map.
11.1002 + template <typename CutMap>
11.1003 + void minCutMap(CutMap& cutMap) const {
11.1004 + for (NodeIt n(_graph); n != INVALID; ++n) {
11.1005 + cutMap.set(n, minCut(n));
11.1006 + }
11.1007 + }
11.1008 +
11.1009 + /// \brief Returns the flow on the edge.
11.1010 + ///
11.1011 + /// Sets the \c flowMap to the flow on the edges. This method can
11.1012 + /// be called after the second phase of algorithm.
11.1013 + Value flow(const Edge& edge) const {
11.1014 + return (*_flow)[edge];
11.1015 + }
11.1016 +
11.1017 + /// @}
11.1018 +
11.1019 + };
11.1020 +
11.1021 +} //namespace lemon
11.1022 +
11.1023 +#endif
12.1 --- a/lemon/preflow.h Wed Nov 14 17:53:08 2007 +0000
12.2 +++ b/lemon/preflow.h Sat Nov 17 20:58:11 2007 +0000
12.3 @@ -19,897 +19,897 @@
12.4 #ifndef LEMON_PREFLOW_H
12.5 #define LEMON_PREFLOW_H
12.6
12.7 -#include <vector>
12.8 -#include <queue>
12.9 -
12.10 #include <lemon/error.h>
12.11 #include <lemon/bits/invalid.h>
12.12 #include <lemon/tolerance.h>
12.13 #include <lemon/maps.h>
12.14 #include <lemon/graph_utils.h>
12.15 +#include <lemon/elevator.h>
12.16
12.17 /// \file
12.18 /// \ingroup max_flow
12.19 /// \brief Implementation of the preflow algorithm.
12.20
12.21 -namespace lemon {
12.22 +namespace lemon {
12.23 +
12.24 + /// \brief Default traits class of Preflow class.
12.25 + ///
12.26 + /// Default traits class of Preflow class.
12.27 + /// \param _Graph Graph type.
12.28 + /// \param _CapacityMap Type of capacity map.
12.29 + template <typename _Graph, typename _CapacityMap>
12.30 + struct PreflowDefaultTraits {
12.31
12.32 - ///\ingroup max_flow
12.33 - ///\brief %Preflow algorithms class.
12.34 + /// \brief The graph type the algorithm runs on.
12.35 + typedef _Graph Graph;
12.36 +
12.37 + /// \brief The type of the map that stores the edge capacities.
12.38 + ///
12.39 + /// The type of the map that stores the edge capacities.
12.40 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
12.41 + typedef _CapacityMap CapacityMap;
12.42 +
12.43 + /// \brief The type of the length of the edges.
12.44 + typedef typename CapacityMap::Value Value;
12.45 +
12.46 + /// \brief The map type that stores the flow values.
12.47 + ///
12.48 + /// The map type that stores the flow values.
12.49 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
12.50 + typedef typename Graph::template EdgeMap<Value> FlowMap;
12.51 +
12.52 + /// \brief Instantiates a FlowMap.
12.53 + ///
12.54 + /// This function instantiates a \ref FlowMap.
12.55 + /// \param graph The graph, to which we would like to define the flow map.
12.56 + static FlowMap* createFlowMap(const Graph& graph) {
12.57 + return new FlowMap(graph);
12.58 + }
12.59 +
12.60 + /// \brief The eleavator type used by Preflow algorithm.
12.61 + ///
12.62 + /// The elevator type used by Preflow algorithm.
12.63 + ///
12.64 + /// \sa Elevator
12.65 + /// \sa LinkedElevator
12.66 + typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
12.67 +
12.68 + /// \brief Instantiates an Elevator.
12.69 + ///
12.70 + /// This function instantiates a \ref Elevator.
12.71 + /// \param graph The graph, to which we would like to define the elevator.
12.72 + /// \param max_level The maximum level of the elevator.
12.73 + static Elevator* createElevator(const Graph& graph, int max_level) {
12.74 + return new Elevator(graph, max_level);
12.75 + }
12.76 +
12.77 + /// \brief The tolerance used by the algorithm
12.78 + ///
12.79 + /// The tolerance used by the algorithm to handle inexact computation.
12.80 + typedef Tolerance<Value> Tolerance;
12.81 +
12.82 + };
12.83 +
12.84 +
12.85 + /// \ingroup max_flow
12.86 ///
12.87 - ///This class provides an implementation of the \e preflow \e
12.88 - ///algorithm producing a flow of maximum value in a directed
12.89 - ///graph. The preflow algorithms are the fastest known max flow algorithms.
12.90 - ///The \e source node, the \e target node, the \e
12.91 - ///capacity of the edges and the \e starting \e flow value of the
12.92 - ///edges should be passed to the algorithm through the
12.93 - ///constructor. It is possible to change these quantities using the
12.94 - ///functions \ref source, \ref target, \ref capacityMap and \ref
12.95 - ///flowMap.
12.96 + /// \brief %Preflow algorithms class.
12.97 ///
12.98 - ///After running \ref lemon::Preflow::phase1() "phase1()"
12.99 - ///or \ref lemon::Preflow::run() "run()", the maximal flow
12.100 - ///value can be obtained by calling \ref flowValue(). The minimum
12.101 - ///value cut can be written into a <tt>bool</tt> node map by
12.102 - ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
12.103 - ///the inclusionwise minimum and maximum of the minimum value cuts,
12.104 - ///resp.)
12.105 + /// This class provides an implementation of the Goldberg's \e
12.106 + /// preflow \e algorithm producing a flow of maximum value in a
12.107 + /// directed graph. The preflow algorithms are the fastest known max
12.108 + /// flow algorithms. The current implementation use a mixture of the
12.109 + /// \e "highest label" and the \e "bound decrease" heuristics.
12.110 + /// The worst case time complexity of the algorithm is \f$O(n^3)\f$.
12.111 ///
12.112 - ///\param Graph The directed graph type the algorithm runs on.
12.113 - ///\param Num The number type of the capacities and the flow values.
12.114 - ///\param CapacityMap The capacity map type.
12.115 - ///\param FlowMap The flow map type.
12.116 - ///\param Tol The tolerance type.
12.117 + /// The algorithm consists from two phases. After the first phase
12.118 + /// the maximal flow value and the minimum cut can be obtained. The
12.119 + /// second phase constructs the feasible maximum flow on each edge.
12.120 ///
12.121 - ///\author Jacint Szabo
12.122 - ///\todo Second template parameter is superfluous
12.123 - template <typename Graph, typename Num,
12.124 - typename CapacityMap=typename Graph::template EdgeMap<Num>,
12.125 - typename FlowMap=typename Graph::template EdgeMap<Num>,
12.126 - typename Tol=Tolerance<Num> >
12.127 + /// \param _Graph The directed graph type the algorithm runs on.
12.128 + /// \param _CapacityMap The flow map type.
12.129 + /// \param _Traits Traits class to set various data types used by
12.130 + /// the algorithm. The default traits class is \ref
12.131 + /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
12.132 + /// documentation of a %Preflow traits class.
12.133 + ///
12.134 + ///\author Jacint Szabo and Balazs Dezso
12.135 +#ifdef DOXYGEN
12.136 + template <typename _Graph, typename _CapacityMap, typename _Traits>
12.137 +#else
12.138 + template <typename _Graph,
12.139 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
12.140 + typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
12.141 +#endif
12.142 class Preflow {
12.143 - protected:
12.144 - typedef typename Graph::Node Node;
12.145 - typedef typename Graph::NodeIt NodeIt;
12.146 - typedef typename Graph::EdgeIt EdgeIt;
12.147 - typedef typename Graph::OutEdgeIt OutEdgeIt;
12.148 - typedef typename Graph::InEdgeIt InEdgeIt;
12.149 + public:
12.150 +
12.151 + typedef _Traits Traits;
12.152 + typedef typename Traits::Graph Graph;
12.153 + typedef typename Traits::CapacityMap CapacityMap;
12.154 + typedef typename Traits::Value Value;
12.155
12.156 - typedef typename Graph::template NodeMap<Node> NNMap;
12.157 - typedef typename std::vector<Node> VecNode;
12.158 + typedef typename Traits::FlowMap FlowMap;
12.159 + typedef typename Traits::Elevator Elevator;
12.160 + typedef typename Traits::Tolerance Tolerance;
12.161
12.162 - const Graph* _g;
12.163 - Node _source;
12.164 - Node _target;
12.165 + /// \brief \ref Exception for uninitialized parameters.
12.166 + ///
12.167 + /// This error represents problems in the initialization
12.168 + /// of the parameters of the algorithms.
12.169 + class UninitializedParameter : public lemon::UninitializedParameter {
12.170 + public:
12.171 + virtual const char* what() const throw() {
12.172 + return "lemon::Preflow::UninitializedParameter";
12.173 + }
12.174 + };
12.175 +
12.176 + private:
12.177 +
12.178 + GRAPH_TYPEDEFS(typename Graph);
12.179 +
12.180 + const Graph& _graph;
12.181 const CapacityMap* _capacity;
12.182 +
12.183 + int _node_num;
12.184 +
12.185 + Node _source, _target;
12.186 +
12.187 FlowMap* _flow;
12.188 + bool _local_flow;
12.189
12.190 - Tol _surely;
12.191 -
12.192 - int _node_num; //the number of nodes of G
12.193 -
12.194 - typename Graph::template NodeMap<int> level;
12.195 - typename Graph::template NodeMap<Num> excess;
12.196 + Elevator* _level;
12.197 + bool _local_level;
12.198
12.199 - // constants used for heuristics
12.200 - static const int H0=20;
12.201 - static const int H1=1;
12.202 + typedef typename Graph::template NodeMap<Value> ExcessMap;
12.203 + ExcessMap* _excess;
12.204 +
12.205 + Tolerance _tolerance;
12.206 +
12.207 + bool _phase;
12.208 +
12.209 + void createStructures() {
12.210 + _node_num = countNodes(_graph);
12.211 +
12.212 + if (!_flow) {
12.213 + _flow = Traits::createFlowMap(_graph);
12.214 + _local_flow = true;
12.215 + }
12.216 + if (!_level) {
12.217 + _level = Traits::createElevator(_graph, _node_num);
12.218 + _local_level = true;
12.219 + }
12.220 + if (!_excess) {
12.221 + _excess = new ExcessMap(_graph);
12.222 + }
12.223 + }
12.224 +
12.225 + void destroyStructures() {
12.226 + if (_local_flow) {
12.227 + delete _flow;
12.228 + }
12.229 + if (_local_level) {
12.230 + delete _level;
12.231 + }
12.232 + if (_excess) {
12.233 + delete _excess;
12.234 + }
12.235 + }
12.236
12.237 public:
12.238
12.239 - ///\ref Exception for the case when s=t.
12.240 + typedef Preflow Create;
12.241
12.242 - ///\ref Exception for the case when the source equals the target.
12.243 - class InvalidArgument : public lemon::LogicError {
12.244 - public:
12.245 - virtual const char* what() const throw() {
12.246 - return "lemon::Preflow::InvalidArgument";
12.247 + ///\name Named template parameters
12.248 +
12.249 + ///@{
12.250 +
12.251 + template <typename _FlowMap>
12.252 + struct DefFlowMapTraits : public Traits {
12.253 + typedef _FlowMap FlowMap;
12.254 + static FlowMap *createFlowMap(const Graph&) {
12.255 + throw UninitializedParameter();
12.256 }
12.257 };
12.258 -
12.259 -
12.260 - ///Indicates the property of the starting flow map.
12.261 -
12.262 - ///Indicates the property of the starting flow map.
12.263 +
12.264 + /// \brief \ref named-templ-param "Named parameter" for setting
12.265 + /// FlowMap type
12.266 ///
12.267 - enum FlowEnum{
12.268 - ///indicates an unspecified edge map. \c flow will be
12.269 - ///set to the constant zero flow in the beginning of
12.270 - ///the algorithm in this case.
12.271 - NO_FLOW,
12.272 - ///constant zero flow
12.273 - ZERO_FLOW,
12.274 - ///any flow, i.e. the sum of the in-flows equals to
12.275 - ///the sum of the out-flows in every node except the \c source and
12.276 - ///the \c target.
12.277 - GEN_FLOW,
12.278 - ///any preflow, i.e. the sum of the in-flows is at
12.279 - ///least the sum of the out-flows in every node except the \c source.
12.280 - PRE_FLOW
12.281 + /// \ref named-templ-param "Named parameter" for setting FlowMap
12.282 + /// type
12.283 + template <typename _FlowMap>
12.284 + struct DefFlowMap
12.285 + : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
12.286 + typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
12.287 };
12.288
12.289 - ///Indicates the state of the preflow algorithm.
12.290 + template <typename _Elevator>
12.291 + struct DefElevatorTraits : public Traits {
12.292 + typedef _Elevator Elevator;
12.293 + static Elevator *createElevator(const Graph&, int) {
12.294 + throw UninitializedParameter();
12.295 + }
12.296 + };
12.297
12.298 - ///Indicates the state of the preflow algorithm.
12.299 + /// \brief \ref named-templ-param "Named parameter" for setting
12.300 + /// Elevator type
12.301 ///
12.302 - enum StatusEnum {
12.303 - ///before running the algorithm or
12.304 - ///at an unspecified state.
12.305 - AFTER_NOTHING,
12.306 - ///right after running \ref phase1()
12.307 - AFTER_PREFLOW_PHASE_1,
12.308 - ///after running \ref phase2()
12.309 - AFTER_PREFLOW_PHASE_2
12.310 + /// \ref named-templ-param "Named parameter" for setting Elevator
12.311 + /// type
12.312 + template <typename _Elevator>
12.313 + struct DefElevator
12.314 + : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
12.315 + typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
12.316 };
12.317 -
12.318 - protected:
12.319 - FlowEnum flow_prop;
12.320 - StatusEnum status; // Do not needle this flag only if necessary.
12.321 -
12.322 - public:
12.323 - ///The constructor of the class.
12.324
12.325 - ///The constructor of the class.
12.326 - ///\param _gr The directed graph the algorithm runs on.
12.327 - ///\param _s The source node.
12.328 - ///\param _t The target node.
12.329 - ///\param _cap The capacity of the edges.
12.330 - ///\param _f The flow of the edges.
12.331 - ///\param _sr Tol class.
12.332 - ///Except the graph, all of these parameters can be reset by
12.333 - ///calling \ref source, \ref target, \ref capacityMap and \ref
12.334 - ///flowMap, resp.
12.335 - Preflow(const Graph& _gr, Node _s, Node _t,
12.336 - const CapacityMap& _cap, FlowMap& _f,
12.337 - const Tol &_sr=Tol()) :
12.338 - _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
12.339 - _flow(&_f), _surely(_sr),
12.340 - _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
12.341 - flow_prop(NO_FLOW), status(AFTER_NOTHING) {
12.342 - if ( _source==_target )
12.343 - throw InvalidArgument();
12.344 - }
12.345 -
12.346 - ///Give a reference to the tolerance handler class
12.347 + template <typename _Elevator>
12.348 + struct DefStandardElevatorTraits : public Traits {
12.349 + typedef _Elevator Elevator;
12.350 + static Elevator *createElevator(const Graph& graph, int max_level) {
12.351 + return new Elevator(graph, max_level);
12.352 + }
12.353 + };
12.354
12.355 - ///Give a reference to the tolerance handler class
12.356 - ///\sa Tolerance
12.357 - Tol &tolerance() { return _surely; }
12.358 + /// \brief \ref named-templ-param "Named parameter" for setting
12.359 + /// Elevator type
12.360 + ///
12.361 + /// \ref named-templ-param "Named parameter" for setting Elevator
12.362 + /// type. The Elevator should be standard constructor interface, ie.
12.363 + /// the graph and the maximum level should be passed to it.
12.364 + template <typename _Elevator>
12.365 + struct DefStandardElevator
12.366 + : public Preflow<Graph, CapacityMap,
12.367 + DefStandardElevatorTraits<_Elevator> > {
12.368 + typedef Preflow<Graph, CapacityMap,
12.369 + DefStandardElevatorTraits<_Elevator> > Create;
12.370 + };
12.371
12.372 - ///Runs the preflow algorithm.
12.373 + /// @}
12.374
12.375 - ///Runs the preflow algorithm.
12.376 + /// \brief The constructor of the class.
12.377 ///
12.378 - void run() {
12.379 - phase1(flow_prop);
12.380 - phase2();
12.381 - }
12.382 -
12.383 - ///Runs the preflow algorithm.
12.384 -
12.385 - ///Runs the preflow algorithm.
12.386 - ///\pre The starting flow map must be
12.387 - /// - a constant zero flow if \c fp is \c ZERO_FLOW,
12.388 - /// - an arbitrary flow if \c fp is \c GEN_FLOW,
12.389 - /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
12.390 - /// - any map if \c fp is NO_FLOW.
12.391 - ///If the starting flow map is a flow or a preflow then
12.392 - ///the algorithm terminates faster.
12.393 - void run(FlowEnum fp) {
12.394 - flow_prop=fp;
12.395 - run();
12.396 - }
12.397 -
12.398 - ///Runs the first phase of the preflow algorithm.
12.399 + /// The constructor of the class.
12.400 + /// \param graph The directed graph the algorithm runs on.
12.401 + /// \param capacity The capacity of the edges.
12.402 + /// \param source The source node.
12.403 + /// \param target The target node.
12.404 + Preflow(const Graph& graph, const CapacityMap& capacity,
12.405 + Node source, Node target)
12.406 + : _graph(graph), _capacity(&capacity),
12.407 + _node_num(0), _source(source), _target(target),
12.408 + _flow(0), _local_flow(false),
12.409 + _level(0), _local_level(false),
12.410 + _excess(0), _tolerance(), _phase() {}
12.411
12.412 - ///The preflow algorithm consists of two phases, this method runs
12.413 - ///the first phase. After the first phase the maximum flow value
12.414 - ///and a minimum value cut can already be computed, although a
12.415 - ///maximum flow is not yet obtained. So after calling this method
12.416 - ///\ref flowValue returns the value of a maximum flow and \ref
12.417 - ///minCut returns a minimum cut.
12.418 - ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
12.419 - ///value cuts unless calling \ref phase2.
12.420 - ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
12.421 - ///is needed for this phase.
12.422 - ///\pre The starting flow must be
12.423 - ///- a constant zero flow if \c fp is \c ZERO_FLOW,
12.424 - ///- an arbitary flow if \c fp is \c GEN_FLOW,
12.425 - ///- an arbitary preflow if \c fp is \c PRE_FLOW,
12.426 - ///- any map if \c fp is NO_FLOW.
12.427 - void phase1(FlowEnum fp)
12.428 - {
12.429 - flow_prop=fp;
12.430 - phase1();
12.431 + /// \brief Destrcutor.
12.432 + ///
12.433 + /// Destructor.
12.434 + ~Preflow() {
12.435 + destroyStructures();
12.436 }
12.437
12.438 + /// \brief Sets the capacity map.
12.439 + ///
12.440 + /// Sets the capacity map.
12.441 + /// \return \c (*this)
12.442 + Preflow& capacityMap(const CapacityMap& map) {
12.443 + _capacity = ↦
12.444 + return *this;
12.445 + }
12.446 +
12.447 + /// \brief Sets the flow map.
12.448 + ///
12.449 + /// Sets the flow map.
12.450 + /// \return \c (*this)
12.451 + Preflow& flowMap(FlowMap& map) {
12.452 + if (_local_flow) {
12.453 + delete _flow;
12.454 + _local_flow = false;
12.455 + }
12.456 + _flow = ↦
12.457 + return *this;
12.458 + }
12.459 +
12.460 + /// \brief Returns the flow map.
12.461 + ///
12.462 + /// \return The flow map.
12.463 + const FlowMap& flowMap() {
12.464 + return *_flow;
12.465 + }
12.466 +
12.467 + /// \brief Sets the elevator.
12.468 + ///
12.469 + /// Sets the elevator.
12.470 + /// \return \c (*this)
12.471 + Preflow& elevator(Elevator& elevator) {
12.472 + if (_local_level) {
12.473 + delete _level;
12.474 + _local_level = false;
12.475 + }
12.476 + _level = &elevator;
12.477 + return *this;
12.478 + }
12.479 +
12.480 + /// \brief Returns the elevator.
12.481 + ///
12.482 + /// \return The elevator.
12.483 + const Elevator& elevator() {
12.484 + return *_level;
12.485 + }
12.486 +
12.487 + /// \brief Sets the source node.
12.488 + ///
12.489 + /// Sets the source node.
12.490 + /// \return \c (*this)
12.491 + Preflow& source(const Node& node) {
12.492 + _source = node;
12.493 + return *this;
12.494 + }
12.495 +
12.496 + /// \brief Sets the target node.
12.497 + ///
12.498 + /// Sets the target node.
12.499 + /// \return \c (*this)
12.500 + Preflow& target(const Node& node) {
12.501 + _target = node;
12.502 + return *this;
12.503 + }
12.504 +
12.505 + /// \brief Sets the tolerance used by algorithm.
12.506 + ///
12.507 + /// Sets the tolerance used by algorithm.
12.508 + Preflow& tolerance(const Tolerance& tolerance) const {
12.509 + _tolerance = tolerance;
12.510 + return *this;
12.511 + }
12.512 +
12.513 + /// \brief Returns the tolerance used by algorithm.
12.514 + ///
12.515 + /// Returns the tolerance used by algorithm.
12.516 + const Tolerance& tolerance() const {
12.517 + return tolerance;
12.518 + }
12.519 +
12.520 + /// \name Execution control The simplest way to execute the
12.521 + /// algorithm is to use one of the member functions called \c
12.522 + /// run().
12.523 + /// \n
12.524 + /// If you need more control on initial solution or
12.525 + /// execution then you have to call one \ref init() function and then
12.526 + /// the startFirstPhase() and if you need the startSecondPhase().
12.527
12.528 - ///Runs the first phase of the preflow algorithm.
12.529 + ///@{
12.530
12.531 - ///The preflow algorithm consists of two phases, this method runs
12.532 - ///the first phase. After the first phase the maximum flow value
12.533 - ///and a minimum value cut can already be computed, although a
12.534 - ///maximum flow is not yet obtained. So after calling this method
12.535 - ///\ref flowValue returns the value of a maximum flow and \ref
12.536 - ///minCut returns a minimum cut.
12.537 - ///\warning \ref minMinCut() and \ref maxMinCut() do not
12.538 - ///give minimum value cuts unless calling \ref phase2().
12.539 - ///\warning A real flow map (i.e. not \ref lemon::NullMap "NullMap")
12.540 - ///is needed for this phase.
12.541 - void phase1()
12.542 - {
12.543 - int heur0=int(H0*_node_num); //time while running 'bound decrease'
12.544 - int heur1=int(H1*_node_num); //time while running 'highest label'
12.545 - int heur=heur1; //starting time interval (#of relabels)
12.546 - int numrelabel=0;
12.547 + /// \brief Initializes the internal data structures.
12.548 + ///
12.549 + /// Initializes the internal data structures.
12.550 + ///
12.551 + void init() {
12.552 + createStructures();
12.553
12.554 - bool what_heur=1;
12.555 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
12.556 + _phase = true;
12.557 + for (NodeIt n(_graph); n != INVALID; ++n) {
12.558 + _excess->set(n, 0);
12.559 + }
12.560
12.561 - bool end=false;
12.562 - //Needed for 'bound decrease', true means no active
12.563 - //nodes are above bound b.
12.564 + for (EdgeIt e(_graph); e != INVALID; ++e) {
12.565 + _flow->set(e, 0);
12.566 + }
12.567
12.568 - int k=_node_num-2; //bound on the highest level under n containing a node
12.569 - int b=k; //bound on the highest level under n containing an active node
12.570 + typename Graph::template NodeMap<bool> reached(_graph, false);
12.571
12.572 - VecNode first(_node_num, INVALID);
12.573 - NNMap next(*_g, INVALID);
12.574 + _level->initStart();
12.575 + _level->initAddItem(_target);
12.576
12.577 - NNMap left(*_g, INVALID);
12.578 - NNMap right(*_g, INVALID);
12.579 - VecNode level_list(_node_num,INVALID);
12.580 - //List of the nodes in level i<n, set to n.
12.581 + std::vector<Node> queue;
12.582 + reached.set(_source, true);
12.583
12.584 - preflowPreproc(first, next, level_list, left, right);
12.585 -
12.586 - //Push/relabel on the highest level active nodes.
12.587 - while ( true ) {
12.588 - if ( b == 0 ) {
12.589 - if ( !what_heur && !end && k > 0 ) {
12.590 - b=k;
12.591 - end=true;
12.592 - } else break;
12.593 - }
12.594 -
12.595 - if ( first[b]==INVALID ) --b;
12.596 - else {
12.597 - end=false;
12.598 - Node w=first[b];
12.599 - first[b]=next[w];
12.600 - int newlevel=push(w, next, first);
12.601 - if ( excess[w] != 0 ) {
12.602 - relabel(w, newlevel, first, next, level_list,
12.603 - left, right, b, k, what_heur);
12.604 - }
12.605 -
12.606 - ++numrelabel;
12.607 - if ( numrelabel >= heur ) {
12.608 - numrelabel=0;
12.609 - if ( what_heur ) {
12.610 - what_heur=0;
12.611 - heur=heur0;
12.612 - end=false;
12.613 - } else {
12.614 - what_heur=1;
12.615 - heur=heur1;
12.616 - b=k;
12.617 + queue.push_back(_target);
12.618 + reached.set(_target, true);
12.619 + while (!queue.empty()) {
12.620 + std::vector<Node> nqueue;
12.621 + for (int i = 0; i < int(queue.size()); ++i) {
12.622 + Node n = queue[i];
12.623 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.624 + Node u = _graph.source(e);
12.625 + if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
12.626 + reached.set(u, true);
12.627 + _level->initAddItem(u);
12.628 + nqueue.push_back(u);
12.629 }
12.630 }
12.631 }
12.632 + queue.swap(nqueue);
12.633 }
12.634 - flow_prop=PRE_FLOW;
12.635 - status=AFTER_PREFLOW_PHASE_1;
12.636 - }
12.637 - // Heuristics:
12.638 - // 2 phase
12.639 - // gap
12.640 - // list 'level_list' on the nodes on level i implemented by hand
12.641 - // stack 'active' on the active nodes on level i
12.642 - // runs heuristic 'highest label' for H1*n relabels
12.643 - // runs heuristic 'bound decrease' for H0*n relabels,
12.644 - // starts with 'highest label'
12.645 - // Parameters H0 and H1 are initialized to 20 and 1.
12.646 + _level->initFinish();
12.647
12.648 -
12.649 - ///Runs the second phase of the preflow algorithm.
12.650 -
12.651 - ///The preflow algorithm consists of two phases, this method runs
12.652 - ///the second phase. After calling \ref phase1() and then
12.653 - ///\ref phase2(),
12.654 - /// \ref flowMap() return a maximum flow, \ref flowValue
12.655 - ///returns the value of a maximum flow, \ref minCut returns a
12.656 - ///minimum cut, while the methods \ref minMinCut and \ref
12.657 - ///maxMinCut return the inclusionwise minimum and maximum cuts of
12.658 - ///minimum value, resp. \pre \ref phase1 must be called before.
12.659 - ///
12.660 - /// \todo The inexact computation can cause positive excess on a set of
12.661 - /// unpushable nodes. We may have to watch the empty level in this case
12.662 - /// due to avoid the terrible long running time.
12.663 - void phase2()
12.664 - {
12.665 -
12.666 - int k=_node_num-2; //bound on the highest level under n containing a node
12.667 - int b=k; //bound on the highest level under n of an active node
12.668 -
12.669 -
12.670 - VecNode first(_node_num, INVALID);
12.671 - NNMap next(*_g, INVALID);
12.672 - level.set(_source,0);
12.673 - std::queue<Node> bfs_queue;
12.674 - bfs_queue.push(_source);
12.675 -
12.676 - while ( !bfs_queue.empty() ) {
12.677 -
12.678 - Node v=bfs_queue.front();
12.679 - bfs_queue.pop();
12.680 - int l=level[v]+1;
12.681 -
12.682 - for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
12.683 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
12.684 - Node u=_g->source(e);
12.685 - if ( level[u] >= _node_num ) {
12.686 - bfs_queue.push(u);
12.687 - level.set(u, l);
12.688 - if ( excess[u] != 0 ) {
12.689 - next.set(u,first[l]);
12.690 - first[l]=u;
12.691 - }
12.692 - }
12.693 - }
12.694 -
12.695 - for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
12.696 - if ( !_surely.positive((*_flow)[e]) ) continue;
12.697 - Node u=_g->target(e);
12.698 - if ( level[u] >= _node_num ) {
12.699 - bfs_queue.push(u);
12.700 - level.set(u, l);
12.701 - if ( excess[u] != 0 ) {
12.702 - next.set(u,first[l]);
12.703 - first[l]=u;
12.704 - }
12.705 - }
12.706 - }
12.707 - }
12.708 - b=_node_num-2;
12.709 -
12.710 - while ( true ) {
12.711 -
12.712 - if ( b == 0 ) break;
12.713 - if ( first[b]==INVALID ) --b;
12.714 - else {
12.715 - Node w=first[b];
12.716 - first[b]=next[w];
12.717 - int newlevel=push(w,next, first);
12.718 -
12.719 - if ( newlevel == _node_num) {
12.720 - excess.set(w, 0);
12.721 - level.set(w,_node_num);
12.722 - }
12.723 - //relabel
12.724 - if ( excess[w] != 0 ) {
12.725 - level.set(w,++newlevel);
12.726 - next.set(w,first[newlevel]);
12.727 - first[newlevel]=w;
12.728 - b=newlevel;
12.729 - }
12.730 - }
12.731 - } // while(true)
12.732 - flow_prop=GEN_FLOW;
12.733 - status=AFTER_PREFLOW_PHASE_2;
12.734 - }
12.735 -
12.736 - /// Returns the value of the maximum flow.
12.737 -
12.738 - /// Returns the value of the maximum flow by returning the excess
12.739 - /// of the target node \c t. This value equals to the value of
12.740 - /// the maximum flow already after running \ref phase1.
12.741 - Num flowValue() const {
12.742 - return excess[_target];
12.743 - }
12.744 -
12.745 -
12.746 - ///Returns a minimum value cut.
12.747 -
12.748 - ///Sets \c M to the characteristic vector of a minimum value
12.749 - ///cut. This method can be called both after running \ref
12.750 - ///phase1 and \ref phase2. It is much faster after
12.751 - ///\ref phase1. \pre M should be a bool-valued node-map. \pre
12.752 - ///If \ref minCut() is called after \ref phase2() then M should
12.753 - ///be initialized to false.
12.754 - template<typename _CutMap>
12.755 - void minCut(_CutMap& M) const {
12.756 - switch ( status ) {
12.757 - case AFTER_PREFLOW_PHASE_1:
12.758 - for(NodeIt v(*_g); v!=INVALID; ++v) {
12.759 - if (level[v] < _node_num) {
12.760 - M.set(v, false);
12.761 - } else {
12.762 - M.set(v, true);
12.763 - }
12.764 - }
12.765 - break;
12.766 - case AFTER_PREFLOW_PHASE_2:
12.767 - minMinCut(M);
12.768 - break;
12.769 - case AFTER_NOTHING:
12.770 - break;
12.771 - }
12.772 - }
12.773 -
12.774 - ///Returns the inclusionwise minimum of the minimum value cuts.
12.775 -
12.776 - ///Sets \c M to the characteristic vector of the minimum value cut
12.777 - ///which is inclusionwise minimum. It is computed by processing a
12.778 - ///bfs from the source node \c s in the residual graph. \pre M
12.779 - ///should be a node map of bools initialized to false. \pre \ref
12.780 - ///phase2 should already be run.
12.781 - template<typename _CutMap>
12.782 - void minMinCut(_CutMap& M) const {
12.783 -
12.784 - std::queue<Node> queue;
12.785 - M.set(_source,true);
12.786 - queue.push(_source);
12.787 -
12.788 - while (!queue.empty()) {
12.789 - Node w=queue.front();
12.790 - queue.pop();
12.791 -
12.792 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.793 - Node v=_g->target(e);
12.794 - if (!M[v] && _surely.positive((*_capacity)[e] -(*_flow)[e]) ) {
12.795 - queue.push(v);
12.796 - M.set(v, true);
12.797 - }
12.798 - }
12.799 -
12.800 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.801 - Node v=_g->source(e);
12.802 - if (!M[v] && _surely.positive((*_flow)[e]) ) {
12.803 - queue.push(v);
12.804 - M.set(v, true);
12.805 - }
12.806 - }
12.807 - }
12.808 - }
12.809 -
12.810 - ///Returns the inclusionwise maximum of the minimum value cuts.
12.811 -
12.812 - ///Sets \c M to the characteristic vector of the minimum value cut
12.813 - ///which is inclusionwise maximum. It is computed by processing a
12.814 - ///backward bfs from the target node \c t in the residual graph.
12.815 - ///\pre \ref phase2() or run() should already be run.
12.816 - template<typename _CutMap>
12.817 - void maxMinCut(_CutMap& M) const {
12.818 -
12.819 - for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
12.820 -
12.821 - std::queue<Node> queue;
12.822 -
12.823 - M.set(_target,false);
12.824 - queue.push(_target);
12.825 -
12.826 - while (!queue.empty()) {
12.827 - Node w=queue.front();
12.828 - queue.pop();
12.829 -
12.830 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.831 - Node v=_g->source(e);
12.832 - if (M[v] && _surely.positive((*_capacity)[e] - (*_flow)[e]) ) {
12.833 - queue.push(v);
12.834 - M.set(v, false);
12.835 - }
12.836 - }
12.837 -
12.838 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.839 - Node v=_g->target(e);
12.840 - if (M[v] && _surely.positive((*_flow)[e]) ) {
12.841 - queue.push(v);
12.842 - M.set(v, false);
12.843 + for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
12.844 + if (_tolerance.positive((*_capacity)[e])) {
12.845 + Node u = _graph.target(e);
12.846 + if ((*_level)[u] == _level->maxLevel()) continue;
12.847 + _flow->set(e, (*_capacity)[e]);
12.848 + _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
12.849 + if (u != _target && !_level->active(u)) {
12.850 + _level->activate(u);
12.851 }
12.852 }
12.853 }
12.854 }
12.855
12.856 - ///Sets the source node to \c _s.
12.857 + /// \brief Initializes the internal data structures.
12.858 + ///
12.859 + /// Initializes the internal data structures and sets the initial
12.860 + /// flow to the given \c flowMap. The \c flowMap should contain a
12.861 + /// flow or at least a preflow, ie. in each node excluding the
12.862 + /// target the incoming flow should greater or equal to the
12.863 + /// outgoing flow.
12.864 + /// \return %False when the given \c flowMap is not a preflow.
12.865 + template <typename FlowMap>
12.866 + bool flowInit(const FlowMap& flowMap) {
12.867 + createStructures();
12.868
12.869 - ///Sets the source node to \c _s.
12.870 - ///
12.871 - void source(Node _s) {
12.872 - _source=_s;
12.873 - if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
12.874 - status=AFTER_NOTHING;
12.875 - }
12.876 + for (EdgeIt e(_graph); e != INVALID; ++e) {
12.877 + _flow->set(e, flowMap[e]);
12.878 + }
12.879
12.880 - ///Returns the source node.
12.881 + for (NodeIt n(_graph); n != INVALID; ++n) {
12.882 + Value excess = 0;
12.883 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.884 + excess += (*_flow)[e];
12.885 + }
12.886 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.887 + excess -= (*_flow)[e];
12.888 + }
12.889 + if (excess < 0 && n != _source) return false;
12.890 + _excess->set(n, excess);
12.891 + }
12.892
12.893 - ///Returns the source node.
12.894 - ///
12.895 - Node source() const {
12.896 - return _source;
12.897 - }
12.898 + typename Graph::template NodeMap<bool> reached(_graph, false);
12.899
12.900 - ///Sets the target node to \c _t.
12.901 + _level->initStart();
12.902 + _level->initAddItem(_target);
12.903
12.904 - ///Sets the target node to \c _t.
12.905 - ///
12.906 - void target(Node _t) {
12.907 - _target=_t;
12.908 - if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
12.909 - status=AFTER_NOTHING;
12.910 - }
12.911 + std::vector<Node> queue;
12.912 + reached.set(_source, true);
12.913
12.914 - ///Returns the target node.
12.915 -
12.916 - ///Returns the target node.
12.917 - ///
12.918 - Node target() const {
12.919 - return _target;
12.920 - }
12.921 -
12.922 - /// Sets the edge map of the capacities to _cap.
12.923 -
12.924 - /// Sets the edge map of the capacities to _cap.
12.925 - ///
12.926 - void capacityMap(const CapacityMap& _cap) {
12.927 - _capacity=&_cap;
12.928 - status=AFTER_NOTHING;
12.929 - }
12.930 - /// Returns a reference to capacity map.
12.931 -
12.932 - /// Returns a reference to capacity map.
12.933 - ///
12.934 - const CapacityMap &capacityMap() const {
12.935 - return *_capacity;
12.936 - }
12.937 -
12.938 - /// Sets the edge map of the flows to _flow.
12.939 -
12.940 - /// Sets the edge map of the flows to _flow.
12.941 - ///
12.942 - void flowMap(FlowMap& _f) {
12.943 - _flow=&_f;
12.944 - flow_prop=NO_FLOW;
12.945 - status=AFTER_NOTHING;
12.946 - }
12.947 -
12.948 - /// Returns a reference to flow map.
12.949 -
12.950 - /// Returns a reference to flow map.
12.951 - ///
12.952 - const FlowMap &flowMap() const {
12.953 - return *_flow;
12.954 - }
12.955 -
12.956 - private:
12.957 -
12.958 - int push(Node w, NNMap& next, VecNode& first) {
12.959 -
12.960 - int lev=level[w];
12.961 - Num exc=excess[w];
12.962 - int newlevel=_node_num; //bound on the next level of w
12.963 -
12.964 - for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.965 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e])) continue;
12.966 - Node v=_g->target(e);
12.967 -
12.968 - if( lev > level[v] ) { //Push is allowed now
12.969 -
12.970 - if ( excess[v] == 0 && v!=_target && v!=_source ) {
12.971 - next.set(v,first[level[v]]);
12.972 - first[level[v]]=v;
12.973 - }
12.974 -
12.975 - Num cap=(*_capacity)[e];
12.976 - Num flo=(*_flow)[e];
12.977 - Num remcap=cap-flo;
12.978 -
12.979 - if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
12.980 -
12.981 - _flow->set(e, flo+exc);
12.982 - excess.set(v, excess[v]+exc);
12.983 - exc=0;
12.984 - break;
12.985 -
12.986 - } else { //A saturating push.
12.987 - _flow->set(e, cap);
12.988 - excess.set(v, excess[v]+remcap);
12.989 - exc-=remcap;
12.990 - }
12.991 - } else if ( newlevel > level[v] ) newlevel = level[v];
12.992 - } //for out edges wv
12.993 -
12.994 - if ( exc != 0 ) {
12.995 - for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
12.996 -
12.997 - if ( !_surely.positive((*_flow)[e]) ) continue;
12.998 - Node v=_g->source(e);
12.999 -
12.1000 - if( lev > level[v] ) { //Push is allowed now
12.1001 -
12.1002 - if ( excess[v] == 0 && v!=_target && v!=_source ) {
12.1003 - next.set(v,first[level[v]]);
12.1004 - first[level[v]]=v;
12.1005 - }
12.1006 -
12.1007 - Num flo=(*_flow)[e];
12.1008 -
12.1009 - if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
12.1010 -
12.1011 - _flow->set(e, flo-exc);
12.1012 - excess.set(v, excess[v]+exc);
12.1013 - exc=0;
12.1014 - break;
12.1015 - } else { //A saturating push.
12.1016 -
12.1017 - excess.set(v, excess[v]+flo);
12.1018 - exc-=flo;
12.1019 - _flow->set(e,0);
12.1020 - }
12.1021 - } else if ( newlevel > level[v] ) newlevel = level[v];
12.1022 - } //for in edges vw
12.1023 -
12.1024 - } // if w still has excess after the out edge for cycle
12.1025 -
12.1026 - excess.set(w, exc);
12.1027 -
12.1028 - return newlevel;
12.1029 - }
12.1030 -
12.1031 -
12.1032 -
12.1033 - void preflowPreproc(VecNode& first, NNMap& next,
12.1034 - VecNode& level_list, NNMap& left, NNMap& right)
12.1035 - {
12.1036 - for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
12.1037 - std::queue<Node> bfs_queue;
12.1038 -
12.1039 - if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
12.1040 - //Reverse_bfs from t in the residual graph,
12.1041 - //to find the starting level.
12.1042 - level.set(_target,0);
12.1043 - bfs_queue.push(_target);
12.1044 -
12.1045 - while ( !bfs_queue.empty() ) {
12.1046 -
12.1047 - Node v=bfs_queue.front();
12.1048 - bfs_queue.pop();
12.1049 - int l=level[v]+1;
12.1050 -
12.1051 - for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
12.1052 - if ( !_surely.positive((*_capacity)[e] - (*_flow)[e] )) continue;
12.1053 - Node w=_g->source(e);
12.1054 - if ( level[w] == _node_num && w != _source ) {
12.1055 - bfs_queue.push(w);
12.1056 - Node z=level_list[l];
12.1057 - if ( z!=INVALID ) left.set(z,w);
12.1058 - right.set(w,z);
12.1059 - level_list[l]=w;
12.1060 - level.set(w, l);
12.1061 + queue.push_back(_target);
12.1062 + reached.set(_target, true);
12.1063 + while (!queue.empty()) {
12.1064 + _level->initNewLevel();
12.1065 + std::vector<Node> nqueue;
12.1066 + for (int i = 0; i < int(queue.size()); ++i) {
12.1067 + Node n = queue[i];
12.1068 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1069 + Node u = _graph.source(e);
12.1070 + if (!reached[u] &&
12.1071 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
12.1072 + reached.set(u, true);
12.1073 + _level->initAddItem(u);
12.1074 + nqueue.push_back(u);
12.1075 }
12.1076 }
12.1077 -
12.1078 - for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
12.1079 - if ( !_surely.positive((*_flow)[e]) ) continue;
12.1080 - Node w=_g->target(e);
12.1081 - if ( level[w] == _node_num && w != _source ) {
12.1082 - bfs_queue.push(w);
12.1083 - Node z=level_list[l];
12.1084 - if ( z!=INVALID ) left.set(z,w);
12.1085 - right.set(w,z);
12.1086 - level_list[l]=w;
12.1087 - level.set(w, l);
12.1088 - }
12.1089 - }
12.1090 - } //while
12.1091 - } //if
12.1092 -
12.1093 -
12.1094 - switch (flow_prop) {
12.1095 - case NO_FLOW:
12.1096 - for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
12.1097 - case ZERO_FLOW:
12.1098 - for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
12.1099 -
12.1100 - //Reverse_bfs from t, to find the starting level.
12.1101 - level.set(_target,0);
12.1102 - bfs_queue.push(_target);
12.1103 -
12.1104 - while ( !bfs_queue.empty() ) {
12.1105 -
12.1106 - Node v=bfs_queue.front();
12.1107 - bfs_queue.pop();
12.1108 - int l=level[v]+1;
12.1109 -
12.1110 - for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
12.1111 - Node w=_g->source(e);
12.1112 - if ( level[w] == _node_num && w != _source ) {
12.1113 - bfs_queue.push(w);
12.1114 - Node z=level_list[l];
12.1115 - if ( z!=INVALID ) left.set(z,w);
12.1116 - right.set(w,z);
12.1117 - level_list[l]=w;
12.1118 - level.set(w, l);
12.1119 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1120 + Node v = _graph.target(e);
12.1121 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
12.1122 + reached.set(v, true);
12.1123 + _level->initAddItem(v);
12.1124 + nqueue.push_back(v);
12.1125 }
12.1126 }
12.1127 }
12.1128 -
12.1129 - //the starting flow
12.1130 - for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
12.1131 - Num c=(*_capacity)[e];
12.1132 - if ( !_surely.positive(c) ) continue;
12.1133 - Node w=_g->target(e);
12.1134 - if ( level[w] < _node_num ) {
12.1135 - if ( excess[w] == 0 && w!=_target ) { //putting into the stack
12.1136 - next.set(w,first[level[w]]);
12.1137 - first[level[w]]=w;
12.1138 - }
12.1139 - _flow->set(e, c);
12.1140 - excess.set(w, excess[w]+c);
12.1141 + queue.swap(nqueue);
12.1142 + }
12.1143 + _level->initFinish();
12.1144 +
12.1145 + for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
12.1146 + Value rem = (*_capacity)[e] - (*_flow)[e];
12.1147 + if (_tolerance.positive(rem)) {
12.1148 + Node u = _graph.target(e);
12.1149 + if ((*_level)[u] == _level->maxLevel()) continue;
12.1150 + _flow->set(e, (*_capacity)[e]);
12.1151 + _excess->set(u, (*_excess)[u] + rem);
12.1152 + if (u != _target && !_level->active(u)) {
12.1153 + _level->activate(u);
12.1154 }
12.1155 }
12.1156 - break;
12.1157 + }
12.1158 + for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
12.1159 + Value rem = (*_flow)[e];
12.1160 + if (_tolerance.positive(rem)) {
12.1161 + Node v = _graph.source(e);
12.1162 + if ((*_level)[v] == _level->maxLevel()) continue;
12.1163 + _flow->set(e, 0);
12.1164 + _excess->set(v, (*_excess)[v] + rem);
12.1165 + if (v != _target && !_level->active(v)) {
12.1166 + _level->activate(v);
12.1167 + }
12.1168 + }
12.1169 + }
12.1170 + return true;
12.1171 + }
12.1172
12.1173 - case GEN_FLOW:
12.1174 - for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
12.1175 - {
12.1176 - Num exc=0;
12.1177 - for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
12.1178 - for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
12.1179 - if (!_surely.positive(exc)) {
12.1180 - exc = 0;
12.1181 - }
12.1182 - excess.set(_target,exc);
12.1183 + /// \brief Starts the first phase of the preflow algorithm.
12.1184 + ///
12.1185 + /// The preflow algorithm consists of two phases, this method runs
12.1186 + /// the first phase. After the first phase the maximum flow value
12.1187 + /// and a minimum value cut can already be computed, although a
12.1188 + /// maximum flow is not yet obtained. So after calling this method
12.1189 + /// \ref flowValue() returns the value of a maximum flow and \ref
12.1190 + /// minCut() returns a minimum cut.
12.1191 + /// \pre One of the \ref init() functions should be called.
12.1192 + void startFirstPhase() {
12.1193 + _phase = true;
12.1194 +
12.1195 + Node n = _level->highestActive();
12.1196 + int level = _level->highestActiveLevel();
12.1197 + while (n != INVALID) {
12.1198 + int num = _node_num;
12.1199 +
12.1200 + while (num > 0 && n != INVALID) {
12.1201 + Value excess = (*_excess)[n];
12.1202 + int new_level = _level->maxLevel();
12.1203 +
12.1204 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1205 + Value rem = (*_capacity)[e] - (*_flow)[e];
12.1206 + if (!_tolerance.positive(rem)) continue;
12.1207 + Node v = _graph.target(e);
12.1208 + if ((*_level)[v] < level) {
12.1209 + if (!_level->active(v) && v != _target) {
12.1210 + _level->activate(v);
12.1211 + }
12.1212 + if (!_tolerance.less(rem, excess)) {
12.1213 + _flow->set(e, (*_flow)[e] + excess);
12.1214 + _excess->set(v, (*_excess)[v] + excess);
12.1215 + excess = 0;
12.1216 + goto no_more_push_1;
12.1217 + } else {
12.1218 + excess -= rem;
12.1219 + _excess->set(v, (*_excess)[v] + rem);
12.1220 + _flow->set(e, (*_capacity)[e]);
12.1221 + }
12.1222 + } else if (new_level > (*_level)[v]) {
12.1223 + new_level = (*_level)[v];
12.1224 + }
12.1225 + }
12.1226 +
12.1227 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1228 + Value rem = (*_flow)[e];
12.1229 + if (!_tolerance.positive(rem)) continue;
12.1230 + Node v = _graph.source(e);
12.1231 + if ((*_level)[v] < level) {
12.1232 + if (!_level->active(v) && v != _target) {
12.1233 + _level->activate(v);
12.1234 + }
12.1235 + if (!_tolerance.less(rem, excess)) {
12.1236 + _flow->set(e, (*_flow)[e] - excess);
12.1237 + _excess->set(v, (*_excess)[v] + excess);
12.1238 + excess = 0;
12.1239 + goto no_more_push_1;
12.1240 + } else {
12.1241 + excess -= rem;
12.1242 + _excess->set(v, (*_excess)[v] + rem);
12.1243 + _flow->set(e, 0);
12.1244 + }
12.1245 + } else if (new_level > (*_level)[v]) {
12.1246 + new_level = (*_level)[v];
12.1247 + }
12.1248 + }
12.1249 +
12.1250 + no_more_push_1:
12.1251 +
12.1252 + _excess->set(n, excess);
12.1253 +
12.1254 + if (excess != 0) {
12.1255 + if (new_level + 1 < _level->maxLevel()) {
12.1256 + _level->liftHighestActive(new_level + 1);
12.1257 + } else {
12.1258 + _level->liftHighestActiveToTop();
12.1259 + }
12.1260 + if (_level->emptyLevel(level)) {
12.1261 + _level->liftToTop(level);
12.1262 + }
12.1263 + } else {
12.1264 + _level->deactivate(n);
12.1265 + }
12.1266 +
12.1267 + n = _level->highestActive();
12.1268 + level = _level->highestActiveLevel();
12.1269 + --num;
12.1270 + }
12.1271 +
12.1272 + num = _node_num * 20;
12.1273 + while (num > 0 && n != INVALID) {
12.1274 + Value excess = (*_excess)[n];
12.1275 + int new_level = _level->maxLevel();
12.1276 +
12.1277 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1278 + Value rem = (*_capacity)[e] - (*_flow)[e];
12.1279 + if (!_tolerance.positive(rem)) continue;
12.1280 + Node v = _graph.target(e);
12.1281 + if ((*_level)[v] < level) {
12.1282 + if (!_level->active(v) && v != _target) {
12.1283 + _level->activate(v);
12.1284 + }
12.1285 + if (!_tolerance.less(rem, excess)) {
12.1286 + _flow->set(e, (*_flow)[e] + excess);
12.1287 + _excess->set(v, (*_excess)[v] + excess);
12.1288 + excess = 0;
12.1289 + goto no_more_push_2;
12.1290 + } else {
12.1291 + excess -= rem;
12.1292 + _excess->set(v, (*_excess)[v] + rem);
12.1293 + _flow->set(e, (*_capacity)[e]);
12.1294 + }
12.1295 + } else if (new_level > (*_level)[v]) {
12.1296 + new_level = (*_level)[v];
12.1297 + }
12.1298 + }
12.1299 +
12.1300 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1301 + Value rem = (*_flow)[e];
12.1302 + if (!_tolerance.positive(rem)) continue;
12.1303 + Node v = _graph.source(e);
12.1304 + if ((*_level)[v] < level) {
12.1305 + if (!_level->active(v) && v != _target) {
12.1306 + _level->activate(v);
12.1307 + }
12.1308 + if (!_tolerance.less(rem, excess)) {
12.1309 + _flow->set(e, (*_flow)[e] - excess);
12.1310 + _excess->set(v, (*_excess)[v] + excess);
12.1311 + excess = 0;
12.1312 + goto no_more_push_2;
12.1313 + } else {
12.1314 + excess -= rem;
12.1315 + _excess->set(v, (*_excess)[v] + rem);
12.1316 + _flow->set(e, 0);
12.1317 + }
12.1318 + } else if (new_level > (*_level)[v]) {
12.1319 + new_level = (*_level)[v];
12.1320 + }
12.1321 + }
12.1322 +
12.1323 + no_more_push_2:
12.1324 +
12.1325 + _excess->set(n, excess);
12.1326 +
12.1327 + if (excess != 0) {
12.1328 + if (new_level + 1 < _level->maxLevel()) {
12.1329 + _level->liftActiveOn(level, new_level + 1);
12.1330 + } else {
12.1331 + _level->liftActiveToTop(level);
12.1332 + }
12.1333 + if (_level->emptyLevel(level)) {
12.1334 + _level->liftToTop(level);
12.1335 + }
12.1336 + } else {
12.1337 + _level->deactivate(n);
12.1338 + }
12.1339 +
12.1340 + while (level >= 0 && _level->activeFree(level)) {
12.1341 + --level;
12.1342 + }
12.1343 + if (level == -1) {
12.1344 + n = _level->highestActive();
12.1345 + level = _level->highestActiveLevel();
12.1346 + } else {
12.1347 + n = _level->activeOn(level);
12.1348 + }
12.1349 + --num;
12.1350 + }
12.1351 + }
12.1352 + }
12.1353 +
12.1354 + /// \brief Starts the second phase of the preflow algorithm.
12.1355 + ///
12.1356 + /// The preflow algorithm consists of two phases, this method runs
12.1357 + /// the second phase. After calling \ref init() and \ref
12.1358 + /// startFirstPhase() and then \ref startSecondPhase(), \ref
12.1359 + /// flowMap() return a maximum flow, \ref flowValue() returns the
12.1360 + /// value of a maximum flow, \ref minCut() returns a minimum cut
12.1361 + /// \pre The \ref init() and startFirstPhase() functions should be
12.1362 + /// called before.
12.1363 + void startSecondPhase() {
12.1364 + _phase = false;
12.1365 +
12.1366 + typename Graph::template NodeMap<bool> reached(_graph, false);
12.1367 + for (NodeIt n(_graph); n != INVALID; ++n) {
12.1368 + reached.set(n, (*_level)[n] < _level->maxLevel());
12.1369 + }
12.1370 +
12.1371 + _level->initStart();
12.1372 + _level->initAddItem(_source);
12.1373 +
12.1374 + std::vector<Node> queue;
12.1375 + queue.push_back(_source);
12.1376 + reached.set(_source, true);
12.1377 +
12.1378 + while (!queue.empty()) {
12.1379 + _level->initNewLevel();
12.1380 + std::vector<Node> nqueue;
12.1381 + for (int i = 0; i < int(queue.size()); ++i) {
12.1382 + Node n = queue[i];
12.1383 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1384 + Node v = _graph.target(e);
12.1385 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
12.1386 + reached.set(v, true);
12.1387 + _level->initAddItem(v);
12.1388 + nqueue.push_back(v);
12.1389 + }
12.1390 + }
12.1391 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1392 + Node u = _graph.source(e);
12.1393 + if (!reached[u] &&
12.1394 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
12.1395 + reached.set(u, true);
12.1396 + _level->initAddItem(u);
12.1397 + nqueue.push_back(u);
12.1398 + }
12.1399 + }
12.1400 + }
12.1401 + queue.swap(nqueue);
12.1402 + }
12.1403 + _level->initFinish();
12.1404 +
12.1405 + for (NodeIt n(_graph); n != INVALID; ++n) {
12.1406 + if ((*_excess)[n] > 0 && _target != n) {
12.1407 + _level->activate(n);
12.1408 + }
12.1409 + }
12.1410 +
12.1411 + Node n;
12.1412 + while ((n = _level->highestActive()) != INVALID) {
12.1413 + Value excess = (*_excess)[n];
12.1414 + int level = _level->highestActiveLevel();
12.1415 + int new_level = _level->maxLevel();
12.1416 +
12.1417 + for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1418 + Value rem = (*_capacity)[e] - (*_flow)[e];
12.1419 + if (!_tolerance.positive(rem)) continue;
12.1420 + Node v = _graph.target(e);
12.1421 + if ((*_level)[v] < level) {
12.1422 + if (!_level->active(v) && v != _source) {
12.1423 + _level->activate(v);
12.1424 + }
12.1425 + if (!_tolerance.less(rem, excess)) {
12.1426 + _flow->set(e, (*_flow)[e] + excess);
12.1427 + _excess->set(v, (*_excess)[v] + excess);
12.1428 + excess = 0;
12.1429 + goto no_more_push;
12.1430 + } else {
12.1431 + excess -= rem;
12.1432 + _excess->set(v, (*_excess)[v] + rem);
12.1433 + _flow->set(e, (*_capacity)[e]);
12.1434 + }
12.1435 + } else if (new_level > (*_level)[v]) {
12.1436 + new_level = (*_level)[v];
12.1437 + }
12.1438 }
12.1439
12.1440 - //the starting flow
12.1441 - for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
12.1442 - Num rem=(*_capacity)[e]-(*_flow)[e];
12.1443 - if ( !_surely.positive(rem) ) continue;
12.1444 - Node w=_g->target(e);
12.1445 - if ( level[w] < _node_num ) {
12.1446 - if ( excess[w] == 0 && w!=_target ) { //putting into the stack
12.1447 - next.set(w,first[level[w]]);
12.1448 - first[level[w]]=w;
12.1449 - }
12.1450 - _flow->set(e, (*_capacity)[e]);
12.1451 - excess.set(w, excess[w]+rem);
12.1452 + for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
12.1453 + Value rem = (*_flow)[e];
12.1454 + if (!_tolerance.positive(rem)) continue;
12.1455 + Node v = _graph.source(e);
12.1456 + if ((*_level)[v] < level) {
12.1457 + if (!_level->active(v) && v != _source) {
12.1458 + _level->activate(v);
12.1459 + }
12.1460 + if (!_tolerance.less(rem, excess)) {
12.1461 + _flow->set(e, (*_flow)[e] - excess);
12.1462 + _excess->set(v, (*_excess)[v] + excess);
12.1463 + excess = 0;
12.1464 + goto no_more_push;
12.1465 + } else {
12.1466 + excess -= rem;
12.1467 + _excess->set(v, (*_excess)[v] + rem);
12.1468 + _flow->set(e, 0);
12.1469 + }
12.1470 + } else if (new_level > (*_level)[v]) {
12.1471 + new_level = (*_level)[v];
12.1472 }
12.1473 }
12.1474 -
12.1475 - for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
12.1476 - if ( !_surely.positive((*_flow)[e]) ) continue;
12.1477 - Node w=_g->source(e);
12.1478 - if ( level[w] < _node_num ) {
12.1479 - if ( excess[w] == 0 && w!=_target ) {
12.1480 - next.set(w,first[level[w]]);
12.1481 - first[level[w]]=w;
12.1482 - }
12.1483 - excess.set(w, excess[w]+(*_flow)[e]);
12.1484 - _flow->set(e, 0);
12.1485 +
12.1486 + no_more_push:
12.1487 +
12.1488 + _excess->set(n, excess);
12.1489 +
12.1490 + if (excess != 0) {
12.1491 + if (new_level + 1 < _level->maxLevel()) {
12.1492 + _level->liftHighestActive(new_level + 1);
12.1493 + } else {
12.1494 + // Calculation error
12.1495 + _level->liftHighestActiveToTop();
12.1496 }
12.1497 - }
12.1498 - break;
12.1499 -
12.1500 - case PRE_FLOW:
12.1501 - //the starting flow
12.1502 - for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
12.1503 - Num rem=(*_capacity)[e]-(*_flow)[e];
12.1504 - if ( !_surely.positive(rem) ) continue;
12.1505 - Node w=_g->target(e);
12.1506 - if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
12.1507 - }
12.1508 -
12.1509 - for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
12.1510 - if ( !_surely.positive((*_flow)[e]) ) continue;
12.1511 - Node w=_g->source(e);
12.1512 - if ( level[w] < _node_num ) _flow->set(e, 0);
12.1513 - }
12.1514 -
12.1515 - //computing the excess
12.1516 - for(NodeIt w(*_g); w!=INVALID; ++w) {
12.1517 - Num exc=0;
12.1518 - for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
12.1519 - for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
12.1520 - if (!_surely.positive(exc)) {
12.1521 - exc = 0;
12.1522 - }
12.1523 - excess.set(w,exc);
12.1524 -
12.1525 - //putting the active nodes into the stack
12.1526 - int lev=level[w];
12.1527 - if ( exc != 0 && lev < _node_num && Node(w) != _target ) {
12.1528 - next.set(w,first[lev]);
12.1529 - first[lev]=w;
12.1530 - }
12.1531 - }
12.1532 - break;
12.1533 - } //switch
12.1534 - } //preflowPreproc
12.1535 -
12.1536 -
12.1537 - void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
12.1538 - VecNode& level_list, NNMap& left,
12.1539 - NNMap& right, int& b, int& k, bool what_heur )
12.1540 - {
12.1541 -
12.1542 - int lev=level[w];
12.1543 -
12.1544 - Node right_n=right[w];
12.1545 - Node left_n=left[w];
12.1546 -
12.1547 - //unlacing starts
12.1548 - if ( right_n!=INVALID ) {
12.1549 - if ( left_n!=INVALID ) {
12.1550 - right.set(left_n, right_n);
12.1551 - left.set(right_n, left_n);
12.1552 + if (_level->emptyLevel(level)) {
12.1553 + // Calculation error
12.1554 + _level->liftToTop(level);
12.1555 + }
12.1556 } else {
12.1557 - level_list[lev]=right_n;
12.1558 - left.set(right_n, INVALID);
12.1559 - }
12.1560 - } else {
12.1561 - if ( left_n!=INVALID ) {
12.1562 - right.set(left_n, INVALID);
12.1563 - } else {
12.1564 - level_list[lev]=INVALID;
12.1565 - }
12.1566 - }
12.1567 - //unlacing ends
12.1568 -
12.1569 - if ( level_list[lev]==INVALID ) {
12.1570 -
12.1571 - //gapping starts
12.1572 - for (int i=lev; i!=k ; ) {
12.1573 - Node v=level_list[++i];
12.1574 - while ( v!=INVALID ) {
12.1575 - level.set(v,_node_num);
12.1576 - v=right[v];
12.1577 - }
12.1578 - level_list[i]=INVALID;
12.1579 - if ( !what_heur ) first[i]=INVALID;
12.1580 + _level->deactivate(n);
12.1581 }
12.1582
12.1583 - level.set(w,_node_num);
12.1584 - b=lev-1;
12.1585 - k=b;
12.1586 - //gapping ends
12.1587 + }
12.1588 + }
12.1589
12.1590 - } else {
12.1591 + /// \brief Runs the preflow algorithm.
12.1592 + ///
12.1593 + /// Runs the preflow algorithm.
12.1594 + /// \note pf.run() is just a shortcut of the following code.
12.1595 + /// \code
12.1596 + /// pf.init();
12.1597 + /// pf.startFirstPhase();
12.1598 + /// pf.startSecondPhase();
12.1599 + /// \endcode
12.1600 + void run() {
12.1601 + init();
12.1602 + startFirstPhase();
12.1603 + startSecondPhase();
12.1604 + }
12.1605
12.1606 - if ( newlevel == _node_num ) level.set(w,_node_num);
12.1607 - else {
12.1608 - level.set(w,++newlevel);
12.1609 - next.set(w,first[newlevel]);
12.1610 - first[newlevel]=w;
12.1611 - if ( what_heur ) b=newlevel;
12.1612 - if ( k < newlevel ) ++k; //now k=newlevel
12.1613 - Node z=level_list[newlevel];
12.1614 - if ( z!=INVALID ) left.set(z,w);
12.1615 - right.set(w,z);
12.1616 - left.set(w,INVALID);
12.1617 - level_list[newlevel]=w;
12.1618 - }
12.1619 + /// \brief Runs the preflow algorithm to compute the minimum cut.
12.1620 + ///
12.1621 + /// Runs the preflow algorithm to compute the minimum cut.
12.1622 + /// \note pf.runMinCut() is just a shortcut of the following code.
12.1623 + /// \code
12.1624 + /// pf.init();
12.1625 + /// pf.startFirstPhase();
12.1626 + /// \endcode
12.1627 + void runMinCut() {
12.1628 + init();
12.1629 + startFirstPhase();
12.1630 + }
12.1631 +
12.1632 + /// @}
12.1633 +
12.1634 + /// \name Query Functions
12.1635 + /// The result of the %Dijkstra algorithm can be obtained using these
12.1636 + /// functions.\n
12.1637 + /// Before the use of these functions,
12.1638 + /// either run() or start() must be called.
12.1639 +
12.1640 + ///@{
12.1641 +
12.1642 + /// \brief Returns the value of the maximum flow.
12.1643 + ///
12.1644 + /// Returns the value of the maximum flow by returning the excess
12.1645 + /// of the target node \c t. This value equals to the value of
12.1646 + /// the maximum flow already after the first phase.
12.1647 + Value flowValue() const {
12.1648 + return (*_excess)[_target];
12.1649 + }
12.1650 +
12.1651 + /// \brief Returns true when the node is on the source side of minimum cut.
12.1652 + ///
12.1653 + /// Returns true when the node is on the source side of minimum
12.1654 + /// cut. This method can be called both after running \ref
12.1655 + /// startFirstPhase() and \ref startSecondPhase().
12.1656 + bool minCut(const Node& node) const {
12.1657 + return ((*_level)[node] == _level->maxLevel()) == _phase;
12.1658 + }
12.1659 +
12.1660 + /// \brief Returns a minimum value cut.
12.1661 + ///
12.1662 + /// Sets the \c cutMap to the characteristic vector of a minimum value
12.1663 + /// cut. This method can be called both after running \ref
12.1664 + /// startFirstPhase() and \ref startSecondPhase(). The result after second
12.1665 + /// phase could be changed slightly if inexact computation is used.
12.1666 + /// \pre The \c cutMap should be a bool-valued node-map.
12.1667 + template <typename CutMap>
12.1668 + void minCutMap(CutMap& cutMap) const {
12.1669 + for (NodeIt n(_graph); n != INVALID; ++n) {
12.1670 + cutMap.set(n, minCut(n));
12.1671 }
12.1672 - } //relabel
12.1673 + }
12.1674
12.1675 - };
12.1676 + /// \brief Returns the flow on the edge.
12.1677 + ///
12.1678 + /// Sets the \c flowMap to the flow on the edges. This method can
12.1679 + /// be called after the second phase of algorithm.
12.1680 + Value flow(const Edge& edge) const {
12.1681 + return (*_flow)[edge];
12.1682 + }
12.1683 +
12.1684 + /// @}
12.1685 + };
12.1686 +}
12.1687
12.1688 - ///\ingroup max_flow
12.1689 - ///\brief Function type interface for Preflow algorithm.
12.1690 - ///
12.1691 - ///Function type interface for Preflow algorithm.
12.1692 - ///\sa Preflow
12.1693 - template<class GR, class CM, class FM>
12.1694 - Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
12.1695 - typename GR::Node source,
12.1696 - typename GR::Node target,
12.1697 - const CM &cap,
12.1698 - FM &flow
12.1699 - )
12.1700 - {
12.1701 - return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
12.1702 - }
12.1703 -
12.1704 -} //namespace lemon
12.1705 -
12.1706 -#endif //LEMON_PREFLOW_H
12.1707 +#endif
13.1 --- a/test/preflow_test.cc Wed Nov 14 17:53:08 2007 +0000
13.2 +++ b/test/preflow_test.cc Sat Nov 17 20:58:11 2007 +0000
13.3 @@ -28,7 +28,7 @@
13.4
13.5 using namespace lemon;
13.6
13.7 -void check_Preflow()
13.8 +void checkPreflow()
13.9 {
13.10 typedef int VType;
13.11 typedef concepts::Graph Graph;
13.12 @@ -37,35 +37,40 @@
13.13 typedef Graph::Edge Edge;
13.14 typedef concepts::ReadMap<Edge,VType> CapMap;
13.15 typedef concepts::ReadWriteMap<Edge,VType> FlowMap;
13.16 - typedef concepts::ReadWriteMap<Node,bool> CutMap;
13.17 + typedef concepts::WriteMap<Node,bool> CutMap;
13.18
13.19 - typedef Preflow<Graph, int, CapMap, FlowMap> PType;
13.20 -
13.21 Graph g;
13.22 Node n;
13.23 + Edge e;
13.24 CapMap cap;
13.25 FlowMap flow;
13.26 CutMap cut;
13.27
13.28 - PType preflow_test(g,n,n,cap,flow);
13.29 + Preflow<Graph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
13.30
13.31 + preflow_test.capacityMap(cap);
13.32 + flow = preflow_test.flowMap();
13.33 + preflow_test.flowMap(flow);
13.34 + preflow_test.source(n);
13.35 + preflow_test.target(n);
13.36 +
13.37 + preflow_test.init();
13.38 + preflow_test.flowInit(cap);
13.39 + preflow_test.startFirstPhase();
13.40 + preflow_test.startSecondPhase();
13.41 preflow_test.run();
13.42 + preflow_test.runMinCut();
13.43 +
13.44 preflow_test.flowValue();
13.45 - preflow_test.source(n);
13.46 - preflow_test.flowMap(flow);
13.47 + preflow_test.minCut(n);
13.48 + preflow_test.minCutMap(cut);
13.49 + preflow_test.flow(e);
13.50
13.51 - preflow_test.phase1(PType::NO_FLOW);
13.52 - preflow_test.minCut(cut);
13.53 -
13.54 - preflow_test.phase2();
13.55 - preflow_test.target(n);
13.56 - preflow_test.capacityMap(cap);
13.57 - preflow_test.minMinCut(cut);
13.58 - preflow_test.maxMinCut(cut);
13.59 }
13.60
13.61 -int cut_value ( SmartGraph& g, SmartGraph::NodeMap<bool>& cut,
13.62 - SmartGraph::EdgeMap<int>& cap) {
13.63 +int cutValue (const SmartGraph& g,
13.64 + const SmartGraph::NodeMap<bool>& cut,
13.65 + const SmartGraph::EdgeMap<int>& cap) {
13.66
13.67 int c=0;
13.68 for(SmartGraph::EdgeIt e(g); e!=INVALID; ++e) {
13.69 @@ -74,6 +79,29 @@
13.70 return c;
13.71 }
13.72
13.73 +bool checkFlow(const SmartGraph& g,
13.74 + const SmartGraph::EdgeMap<int>& flow,
13.75 + const SmartGraph::EdgeMap<int>& cap,
13.76 + SmartGraph::Node s, SmartGraph::Node t) {
13.77 +
13.78 + for (SmartGraph::EdgeIt e(g); e != INVALID; ++e) {
13.79 + if (flow[e] < 0 || flow[e] > cap[e]) return false;
13.80 + }
13.81 +
13.82 + for (SmartGraph::NodeIt n(g); n != INVALID; ++n) {
13.83 + if (n == s || n == t) continue;
13.84 + int sum = 0;
13.85 + for (SmartGraph::OutEdgeIt e(g, n); e != INVALID; ++e) {
13.86 + sum += flow[e];
13.87 + }
13.88 + for (SmartGraph::InEdgeIt e(g, n); e != INVALID; ++e) {
13.89 + sum -= flow[e];
13.90 + }
13.91 + if (sum != 0) return false;
13.92 + }
13.93 + return true;
13.94 +}
13.95 +
13.96 int main() {
13.97
13.98 typedef SmartGraph Graph;
13.99 @@ -85,7 +113,7 @@
13.100 typedef Graph::EdgeMap<int> FlowMap;
13.101 typedef Graph::NodeMap<bool> CutMap;
13.102
13.103 - typedef Preflow<Graph, int> PType;
13.104 + typedef Preflow<Graph, CapMap> PType;
13.105
13.106 std::string f_name;
13.107 if( getenv("srcdir") )
13.108 @@ -102,74 +130,50 @@
13.109 CapMap cap(g);
13.110 readDimacs(file, g, cap, s, t);
13.111
13.112 - FlowMap flow(g,0);
13.113 + PType preflow_test(g, cap, s, t);
13.114 + preflow_test.run();
13.115 +
13.116 + check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
13.117 + "The flow is not feasible.");
13.118
13.119 -
13.120 + CutMap min_cut(g);
13.121 + preflow_test.minCutMap(min_cut);
13.122 + int min_cut_value=cutValue(g,min_cut,cap);
13.123 +
13.124 + check(preflow_test.flowValue() == min_cut_value,
13.125 + "The max flow value is not equal to the three min cut values.");
13.126
13.127 - PType preflow_test(g, s, t, cap, flow);
13.128 - preflow_test.run(PType::ZERO_FLOW);
13.129 -
13.130 - CutMap min_cut(g,false);
13.131 - preflow_test.minCut(min_cut);
13.132 - int min_cut_value=cut_value(g,min_cut,cap);
13.133 -
13.134 - CutMap min_min_cut(g,false);
13.135 - preflow_test.minMinCut(min_min_cut);
13.136 - int min_min_cut_value=cut_value(g,min_min_cut,cap);
13.137 -
13.138 - CutMap max_min_cut(g,false);
13.139 - preflow_test.maxMinCut(max_min_cut);
13.140 - int max_min_cut_value=cut_value(g,max_min_cut,cap);
13.141 -
13.142 - check(preflow_test.flowValue() == min_cut_value &&
13.143 - min_cut_value == min_min_cut_value &&
13.144 - min_min_cut_value == max_min_cut_value,
13.145 - "The max flow value is not equal to the three min cut values.");
13.146 + FlowMap flow(g);
13.147 + flow = preflow_test.flowMap();
13.148
13.149 int flow_value=preflow_test.flowValue();
13.150
13.151 + for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
13.152 + preflow_test.flowInit(flow);
13.153 + preflow_test.startFirstPhase();
13.154
13.155 -
13.156 - for(EdgeIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
13.157 - preflow_test.capacityMap(cap);
13.158 -
13.159 - preflow_test.phase1(PType::PRE_FLOW);
13.160 -
13.161 - CutMap min_cut1(g,false);
13.162 - preflow_test.minCut(min_cut1);
13.163 - min_cut_value=cut_value(g,min_cut1,cap);
13.164 + CutMap min_cut1(g);
13.165 + preflow_test.minCutMap(min_cut1);
13.166 + min_cut_value=cutValue(g,min_cut1,cap);
13.167
13.168 check(preflow_test.flowValue() == min_cut_value &&
13.169 min_cut_value == 2*flow_value,
13.170 "The max flow value or the min cut value is wrong.");
13.171
13.172 - preflow_test.phase2();
13.173 + preflow_test.startSecondPhase();
13.174
13.175 - CutMap min_cut2(g,false);
13.176 - preflow_test.minCut(min_cut2);
13.177 - min_cut_value=cut_value(g,min_cut2,cap);
13.178 + check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
13.179 + "The flow is not feasible.");
13.180 +
13.181 + CutMap min_cut2(g);
13.182 + preflow_test.minCutMap(min_cut2);
13.183 + min_cut_value=cutValue(g,min_cut2,cap);
13.184
13.185 - CutMap min_min_cut2(g,false);
13.186 - preflow_test.minMinCut(min_min_cut2);
13.187 - min_min_cut_value=cut_value(g,min_min_cut2,cap);
13.188 -
13.189 - preflow_test.maxMinCut(max_min_cut);
13.190 - max_min_cut_value=cut_value(g,max_min_cut,cap);
13.191 -
13.192 check(preflow_test.flowValue() == min_cut_value &&
13.193 - min_cut_value == min_min_cut_value &&
13.194 - min_min_cut_value == max_min_cut_value &&
13.195 min_cut_value == 2*flow_value,
13.196 "The max flow value or the three min cut values were not doubled");
13.197
13.198
13.199 -
13.200 - EdgeIt e(g);
13.201 - for( int i=1; i==10; ++i ) {
13.202 - flow.set(e,0);
13.203 - ++e;
13.204 - }
13.205 -
13.206 preflow_test.flowMap(flow);
13.207
13.208 NodeIt tmp1(g,s);
13.209 @@ -185,19 +189,13 @@
13.210
13.211 preflow_test.run();
13.212
13.213 - CutMap min_cut3(g,false);
13.214 - preflow_test.minCut(min_cut3);
13.215 - min_cut_value=cut_value(g,min_cut3,cap);
13.216 + CutMap min_cut3(g);
13.217 + preflow_test.minCutMap(min_cut3);
13.218 + min_cut_value=cutValue(g,min_cut3,cap);
13.219
13.220 - CutMap min_min_cut3(g,false);
13.221 - preflow_test.minMinCut(min_min_cut3);
13.222 - min_min_cut_value=cut_value(g,min_min_cut3,cap);
13.223 -
13.224 - preflow_test.maxMinCut(max_min_cut);
13.225 - max_min_cut_value=cut_value(g,max_min_cut,cap);
13.226
13.227 - check(preflow_test.flowValue() == min_cut_value &&
13.228 - min_cut_value == min_min_cut_value &&
13.229 - min_min_cut_value == max_min_cut_value,
13.230 + check(preflow_test.flowValue() == min_cut_value,
13.231 "The max flow value or the three min cut values are incorrect.");
13.232 +
13.233 + return 0;
13.234 }
