epoc32/include/stdapis/boost/graph/push_relabel_max_flow.hpp
branchSymbian2
changeset 2 2fe1408b6811
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/epoc32/include/stdapis/boost/graph/push_relabel_max_flow.hpp	Tue Mar 16 16:12:26 2010 +0000
     1.3 @@ -0,0 +1,727 @@
     1.4 +//=======================================================================
     1.5 +// Copyright 2000 University of Notre Dame.
     1.6 +// Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee
     1.7 +//
     1.8 +// Distributed under the Boost Software License, Version 1.0. (See
     1.9 +// accompanying file LICENSE_1_0.txt or copy at
    1.10 +// http://www.boost.org/LICENSE_1_0.txt)
    1.11 +//=======================================================================
    1.12 +
    1.13 +#ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP
    1.14 +#define BOOST_PUSH_RELABEL_MAX_FLOW_HPP
    1.15 +
    1.16 +#include <boost/config.hpp>
    1.17 +#include <cassert>
    1.18 +#include <vector>
    1.19 +#include <list>
    1.20 +#include <iosfwd>
    1.21 +#include <algorithm> // for std::min and std::max
    1.22 +
    1.23 +#include <boost/pending/queue.hpp>
    1.24 +#include <boost/limits.hpp>
    1.25 +#include <boost/graph/graph_concepts.hpp>
    1.26 +#include <boost/graph/named_function_params.hpp>
    1.27 +
    1.28 +namespace boost {
    1.29 +
    1.30 +  namespace detail {
    1.31 +    
    1.32 +   // This implementation is based on Goldberg's 
    1.33 +   // "On Implementing Push-Relabel Method for the Maximum Flow Problem"
    1.34 +   // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171
    1.35 +   // and on the h_prf.c and hi_pr.c code written by the above authors.
    1.36 +
    1.37 +   // This implements the highest-label version of the push-relabel method
    1.38 +   // with the global relabeling and gap relabeling heuristics.
    1.39 +
    1.40 +   // The terms "rank", "distance", "height" are synonyms in
    1.41 +   // Goldberg's implementation, paper and in the CLR.  A "layer" is a
    1.42 +   // group of vertices with the same distance. The vertices in each
    1.43 +   // layer are categorized as active or inactive.  An active vertex
    1.44 +   // has positive excess flow and its distance is less than n (it is
    1.45 +   // not blocked).
    1.46 +
    1.47 +    template <class Vertex>
    1.48 +    struct preflow_layer {
    1.49 +      std::list<Vertex> active_vertices;
    1.50 +      std::list<Vertex> inactive_vertices;
    1.51 +    };
    1.52 +
    1.53 +    template <class Graph, 
    1.54 +              class EdgeCapacityMap,    // integer value type
    1.55 +              class ResidualCapacityEdgeMap,
    1.56 +              class ReverseEdgeMap,
    1.57 +              class VertexIndexMap,     // vertex_descriptor -> integer
    1.58 +              class FlowValue>
    1.59 +    class push_relabel
    1.60 +    {
    1.61 +    public:
    1.62 +      typedef graph_traits<Graph> Traits;
    1.63 +      typedef typename Traits::vertex_descriptor vertex_descriptor;
    1.64 +      typedef typename Traits::edge_descriptor edge_descriptor;
    1.65 +      typedef typename Traits::vertex_iterator vertex_iterator;
    1.66 +      typedef typename Traits::out_edge_iterator out_edge_iterator;
    1.67 +      typedef typename Traits::vertices_size_type vertices_size_type;
    1.68 +      typedef typename Traits::edges_size_type edges_size_type;
    1.69 +
    1.70 +      typedef preflow_layer<vertex_descriptor> Layer;
    1.71 +      typedef std::vector< Layer > LayerArray;
    1.72 +      typedef typename LayerArray::iterator layer_iterator;
    1.73 +      typedef typename LayerArray::size_type distance_size_type;
    1.74 +
    1.75 +      typedef color_traits<default_color_type> ColorTraits;
    1.76 +
    1.77 +      //=======================================================================
    1.78 +      // Some helper predicates
    1.79 +
    1.80 +      inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {
    1.81 +        return distance[u] == distance[v] + 1;
    1.82 +      }
    1.83 +      inline bool is_residual_edge(edge_descriptor a) {
    1.84 +        return 0 < residual_capacity[a];
    1.85 +      }
    1.86 +      inline bool is_saturated(edge_descriptor a) {
    1.87 +        return residual_capacity[a] == 0;
    1.88 +      }
    1.89 +
    1.90 +      //=======================================================================
    1.91 +      // Layer List Management Functions
    1.92 +
    1.93 +      typedef typename std::list<vertex_descriptor>::iterator list_iterator;
    1.94 +
    1.95 +      void add_to_active_list(vertex_descriptor u, Layer& layer) {
    1.96 +        BOOST_USING_STD_MIN();
    1.97 +        BOOST_USING_STD_MAX();
    1.98 +        layer.active_vertices.push_front(u);
    1.99 +        max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], max_active);
   1.100 +        min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active);
   1.101 +        layer_list_ptr[u] = layer.active_vertices.begin();
   1.102 +      }
   1.103 +      void remove_from_active_list(vertex_descriptor u) {
   1.104 +        layers[distance[u]].active_vertices.erase(layer_list_ptr[u]);    
   1.105 +      }
   1.106 +
   1.107 +      void add_to_inactive_list(vertex_descriptor u, Layer& layer) {
   1.108 +        layer.inactive_vertices.push_front(u);
   1.109 +        layer_list_ptr[u] = layer.inactive_vertices.begin();
   1.110 +      }
   1.111 +      void remove_from_inactive_list(vertex_descriptor u) {
   1.112 +        layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]);    
   1.113 +      }
   1.114 +
   1.115 +      //=======================================================================
   1.116 +      // initialization
   1.117 +      push_relabel(Graph& g_, 
   1.118 +                   EdgeCapacityMap cap,
   1.119 +                   ResidualCapacityEdgeMap res,
   1.120 +                   ReverseEdgeMap rev,
   1.121 +                   vertex_descriptor src_, 
   1.122 +                   vertex_descriptor sink_,
   1.123 +                   VertexIndexMap idx)
   1.124 +        : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_), 
   1.125 +          index(idx),
   1.126 +          excess_flow(num_vertices(g_)),
   1.127 +          current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second),
   1.128 +          distance(num_vertices(g_)),
   1.129 +          color(num_vertices(g_)),
   1.130 +          reverse_edge(rev),
   1.131 +          residual_capacity(res),
   1.132 +          layers(num_vertices(g_)),
   1.133 +          layer_list_ptr(num_vertices(g_), 
   1.134 +                         layers.front().inactive_vertices.end()),
   1.135 +          push_count(0), update_count(0), relabel_count(0), 
   1.136 +          gap_count(0), gap_node_count(0),
   1.137 +          work_since_last_update(0)
   1.138 +      {
   1.139 +        vertex_iterator u_iter, u_end;
   1.140 +        // Don't count the reverse edges
   1.141 +        edges_size_type m = num_edges(g) / 2;
   1.142 +        nm = alpha() * n + m;
   1.143 +
   1.144 +        // Initialize flow to zero which means initializing
   1.145 +        // the residual capacity to equal the capacity.
   1.146 +        out_edge_iterator ei, e_end;
   1.147 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
   1.148 +          for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) {
   1.149 +            residual_capacity[*ei] = capacity[*ei];
   1.150 +          }
   1.151 +
   1.152 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.153 +          vertex_descriptor u = *u_iter;
   1.154 +          excess_flow[u] = 0;
   1.155 +          current[u] = out_edges(u, g).first;
   1.156 +        }
   1.157 +
   1.158 +        bool overflow_detected = false;
   1.159 +        FlowValue test_excess = 0;
   1.160 +
   1.161 +        out_edge_iterator a_iter, a_end;
   1.162 +        for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter)
   1.163 +          if (target(*a_iter, g) != src)
   1.164 +            test_excess += residual_capacity[*a_iter];
   1.165 +        if (test_excess > (std::numeric_limits<FlowValue>::max)())
   1.166 +          overflow_detected = true;
   1.167 +
   1.168 +        if (overflow_detected)
   1.169 +          excess_flow[src] = (std::numeric_limits<FlowValue>::max)();
   1.170 +        else {
   1.171 +          excess_flow[src] = 0;
   1.172 +          for (tie(a_iter, a_end) = out_edges(src, g); 
   1.173 +               a_iter != a_end; ++a_iter) {
   1.174 +            edge_descriptor a = *a_iter;
   1.175 +            if (target(a, g) != src) {
   1.176 +              ++push_count;
   1.177 +              FlowValue delta = residual_capacity[a];
   1.178 +              residual_capacity[a] -= delta;
   1.179 +              residual_capacity[reverse_edge[a]] += delta;
   1.180 +              excess_flow[target(a, g)] += delta;
   1.181 +            }
   1.182 +          }
   1.183 +        }
   1.184 +        max_distance = num_vertices(g) - 1;
   1.185 +        max_active = 0;
   1.186 +        min_active = n;
   1.187 +
   1.188 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.189 +          vertex_descriptor u = *u_iter;
   1.190 +          if (u == sink) {
   1.191 +            distance[u] = 0;
   1.192 +            continue;
   1.193 +          } else if (u == src && !overflow_detected)
   1.194 +            distance[u] = n;
   1.195 +          else
   1.196 +            distance[u] = 1;
   1.197 +
   1.198 +          if (excess_flow[u] > 0)
   1.199 +            add_to_active_list(u, layers[1]);
   1.200 +          else if (distance[u] < n)
   1.201 +            add_to_inactive_list(u, layers[1]);
   1.202 +        }       
   1.203 +
   1.204 +      } // push_relabel constructor
   1.205 +
   1.206 +      //=======================================================================
   1.207 +      // This is a breadth-first search over the residual graph
   1.208 +      // (well, actually the reverse of the residual graph).
   1.209 +      // Would be cool to have a graph view adaptor for hiding certain
   1.210 +      // edges, like the saturated (non-residual) edges in this case.
   1.211 +      // Goldberg's implementation abused "distance" for the coloring.
   1.212 +      void global_distance_update()
   1.213 +      {
   1.214 +        BOOST_USING_STD_MAX();
   1.215 +        ++update_count;
   1.216 +        vertex_iterator u_iter, u_end;
   1.217 +        for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.218 +          color[*u_iter] = ColorTraits::white();
   1.219 +          distance[*u_iter] = n;
   1.220 +        }
   1.221 +        color[sink] = ColorTraits::gray();
   1.222 +        distance[sink] = 0;
   1.223 +        
   1.224 +        for (distance_size_type l = 0; l <= max_distance; ++l) {
   1.225 +          layers[l].active_vertices.clear();
   1.226 +          layers[l].inactive_vertices.clear();
   1.227 +        }
   1.228 +        
   1.229 +        max_distance = max_active = 0;
   1.230 +        min_active = n;
   1.231 +
   1.232 +        Q.push(sink);
   1.233 +        while (! Q.empty()) {
   1.234 +          vertex_descriptor u = Q.top();
   1.235 +          Q.pop();
   1.236 +          distance_size_type d_v = distance[u] + 1;
   1.237 +
   1.238 +          out_edge_iterator ai, a_end;
   1.239 +          for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
   1.240 +            edge_descriptor a = *ai;
   1.241 +            vertex_descriptor v = target(a, g);
   1.242 +            if (color[v] == ColorTraits::white()
   1.243 +                && is_residual_edge(reverse_edge[a])) {
   1.244 +              distance[v] = d_v;
   1.245 +              color[v] = ColorTraits::gray();
   1.246 +              current[v] = out_edges(v, g).first;
   1.247 +              max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance);
   1.248 +
   1.249 +              if (excess_flow[v] > 0)
   1.250 +                add_to_active_list(v, layers[d_v]);
   1.251 +              else
   1.252 +                add_to_inactive_list(v, layers[d_v]);
   1.253 +
   1.254 +              Q.push(v);
   1.255 +            }
   1.256 +          }
   1.257 +        }
   1.258 +      } // global_distance_update()
   1.259 +
   1.260 +      //=======================================================================
   1.261 +      // This function is called "push" in Goldberg's h_prf implementation,
   1.262 +      // but it is called "discharge" in the paper and in hi_pr.c.
   1.263 +      void discharge(vertex_descriptor u)
   1.264 +      {
   1.265 +        assert(excess_flow[u] > 0);
   1.266 +        while (1) {
   1.267 +          out_edge_iterator ai, ai_end;
   1.268 +          for (ai = current[u], ai_end = out_edges(u, g).second;
   1.269 +               ai != ai_end; ++ai) {
   1.270 +            edge_descriptor a = *ai;
   1.271 +            if (is_residual_edge(a)) {
   1.272 +              vertex_descriptor v = target(a, g);
   1.273 +              if (is_admissible(u, v)) {
   1.274 +                ++push_count;
   1.275 +                if (v != sink && excess_flow[v] == 0) {
   1.276 +                  remove_from_inactive_list(v);
   1.277 +                  add_to_active_list(v, layers[distance[v]]);
   1.278 +                }
   1.279 +                push_flow(a);
   1.280 +                if (excess_flow[u] == 0)
   1.281 +                  break;
   1.282 +              } 
   1.283 +            } 
   1.284 +          } // for out_edges of i starting from current
   1.285 +
   1.286 +          Layer& layer = layers[distance[u]];
   1.287 +          distance_size_type du = distance[u];
   1.288 +
   1.289 +          if (ai == ai_end) {   // i must be relabeled
   1.290 +            relabel_distance(u);
   1.291 +            if (layer.active_vertices.empty()
   1.292 +                && layer.inactive_vertices.empty())
   1.293 +              gap(du);
   1.294 +            if (distance[u] == n)
   1.295 +              break;
   1.296 +          } else {              // i is no longer active
   1.297 +            current[u] = ai;
   1.298 +            add_to_inactive_list(u, layer);
   1.299 +            break;
   1.300 +          }
   1.301 +        } // while (1)
   1.302 +      } // discharge()
   1.303 +
   1.304 +      //=======================================================================
   1.305 +      // This corresponds to the "push" update operation of the paper,
   1.306 +      // not the "push" function in Goldberg's h_prf.c implementation.
   1.307 +      // The idea is to push the excess flow from from vertex u to v.
   1.308 +      void push_flow(edge_descriptor u_v)
   1.309 +      {
   1.310 +        vertex_descriptor
   1.311 +          u = source(u_v, g),
   1.312 +          v = target(u_v, g);
   1.313 +        
   1.314 +        BOOST_USING_STD_MIN();
   1.315 +        FlowValue flow_delta
   1.316 +          = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]);
   1.317 +
   1.318 +        residual_capacity[u_v] -= flow_delta;
   1.319 +        residual_capacity[reverse_edge[u_v]] += flow_delta;
   1.320 +
   1.321 +        excess_flow[u] -= flow_delta;
   1.322 +        excess_flow[v] += flow_delta;
   1.323 +      } // push_flow()
   1.324 +
   1.325 +      //=======================================================================
   1.326 +      // The main purpose of this routine is to set distance[v]
   1.327 +      // to the smallest value allowed by the valid labeling constraints,
   1.328 +      // which are:
   1.329 +      // distance[t] = 0
   1.330 +      // distance[u] <= distance[v] + 1   for every residual edge (u,v)
   1.331 +      //
   1.332 +      distance_size_type relabel_distance(vertex_descriptor u)
   1.333 +      {
   1.334 +        BOOST_USING_STD_MAX();
   1.335 +        ++relabel_count;
   1.336 +        work_since_last_update += beta();
   1.337 +
   1.338 +        distance_size_type min_distance = num_vertices(g);
   1.339 +        distance[u] = min_distance;
   1.340 +
   1.341 +        // Examine the residual out-edges of vertex i, choosing the
   1.342 +        // edge whose target vertex has the minimal distance.
   1.343 +        out_edge_iterator ai, a_end, min_edge_iter;
   1.344 +        for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
   1.345 +          ++work_since_last_update;
   1.346 +          edge_descriptor a = *ai;
   1.347 +          vertex_descriptor v = target(a, g);
   1.348 +          if (is_residual_edge(a) && distance[v] < min_distance) {
   1.349 +            min_distance = distance[v];
   1.350 +            min_edge_iter = ai;
   1.351 +          }
   1.352 +        }
   1.353 +        ++min_distance;
   1.354 +        if (min_distance < n) {
   1.355 +          distance[u] = min_distance;     // this is the main action
   1.356 +          current[u] = min_edge_iter;
   1.357 +          max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance);
   1.358 +        }
   1.359 +        return min_distance;
   1.360 +      } // relabel_distance()
   1.361 +
   1.362 +      //=======================================================================
   1.363 +      // cleanup beyond the gap
   1.364 +      void gap(distance_size_type empty_distance)
   1.365 +      {
   1.366 +        ++gap_count;
   1.367 +
   1.368 +        distance_size_type r; // distance of layer before the current layer
   1.369 +        r = empty_distance - 1;
   1.370 +
   1.371 +        // Set the distance for the vertices beyond the gap to "infinity".
   1.372 +        for (layer_iterator l = layers.begin() + empty_distance + 1;
   1.373 +             l < layers.begin() + max_distance; ++l) {
   1.374 +          list_iterator i;
   1.375 +          for (i = l->inactive_vertices.begin(); 
   1.376 +               i != l->inactive_vertices.end(); ++i) {
   1.377 +            distance[*i] = n;
   1.378 +            ++gap_node_count;
   1.379 +          }
   1.380 +          l->inactive_vertices.clear();
   1.381 +        }
   1.382 +        max_distance = r;
   1.383 +        max_active = r;
   1.384 +      }
   1.385 +
   1.386 +      //=======================================================================
   1.387 +      // This is the core part of the algorithm, "phase one".
   1.388 +      FlowValue maximum_preflow()
   1.389 +      {
   1.390 +        work_since_last_update = 0;
   1.391 +
   1.392 +        while (max_active >= min_active) { // "main" loop
   1.393 +
   1.394 +          Layer& layer = layers[max_active];
   1.395 +          list_iterator u_iter = layer.active_vertices.begin();
   1.396 +
   1.397 +          if (u_iter == layer.active_vertices.end())
   1.398 +            --max_active;
   1.399 +          else {
   1.400 +            vertex_descriptor u = *u_iter;
   1.401 +            remove_from_active_list(u);
   1.402 +            
   1.403 +            discharge(u);
   1.404 +
   1.405 +            if (work_since_last_update * global_update_frequency() > nm) {
   1.406 +              global_distance_update();
   1.407 +              work_since_last_update = 0;
   1.408 +            }
   1.409 +          }
   1.410 +        } // while (max_active >= min_active)
   1.411 +
   1.412 +        return excess_flow[sink];
   1.413 +      } // maximum_preflow()
   1.414 +
   1.415 +      //=======================================================================
   1.416 +      // remove excess flow, the "second phase"
   1.417 +      // This does a DFS on the reverse flow graph of nodes with excess flow.
   1.418 +      // If a cycle is found, cancel it.
   1.419 +      // Return the nodes with excess flow in topological order.
   1.420 +      //
   1.421 +      // Unlike the prefl_to_flow() implementation, we use
   1.422 +      //   "color" instead of "distance" for the DFS labels
   1.423 +      //   "parent" instead of nl_prev for the DFS tree
   1.424 +      //   "topo_next" instead of nl_next for the topological ordering
   1.425 +      void convert_preflow_to_flow()
   1.426 +      {
   1.427 +        vertex_iterator u_iter, u_end;
   1.428 +        out_edge_iterator ai, a_end;
   1.429 +
   1.430 +        vertex_descriptor r, restart, u;
   1.431 +
   1.432 +        std::vector<vertex_descriptor> parent(n);
   1.433 +        std::vector<vertex_descriptor> topo_next(n);
   1.434 +
   1.435 +        vertex_descriptor tos(parent[0]), 
   1.436 +          bos(parent[0]); // bogus initialization, just to avoid warning
   1.437 +        bool bos_null = true;
   1.438 +
   1.439 +        // handle self-loops
   1.440 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
   1.441 +          for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai)
   1.442 +            if (target(*ai, g) == *u_iter)
   1.443 +              residual_capacity[*ai] = capacity[*ai];
   1.444 +
   1.445 +        // initialize
   1.446 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.447 +          u = *u_iter;
   1.448 +          color[u] = ColorTraits::white();
   1.449 +          parent[u] = u;
   1.450 +          current[u] = out_edges(u, g).first;
   1.451 +        }
   1.452 +        // eliminate flow cycles and topologically order the vertices
   1.453 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.454 +          u = *u_iter;
   1.455 +          if (color[u] == ColorTraits::white() 
   1.456 +              && excess_flow[u] > 0
   1.457 +              && u != src && u != sink ) {
   1.458 +            r = u;
   1.459 +            color[r] = ColorTraits::gray();
   1.460 +            while (1) {
   1.461 +              for (; current[u] != out_edges(u, g).second; ++current[u]) {
   1.462 +                edge_descriptor a = *current[u];
   1.463 +                if (capacity[a] == 0 && is_residual_edge(a)) {
   1.464 +                  vertex_descriptor v = target(a, g);
   1.465 +                  if (color[v] == ColorTraits::white()) {
   1.466 +                    color[v] = ColorTraits::gray();
   1.467 +                    parent[v] = u;
   1.468 +                    u = v;
   1.469 +                    break;
   1.470 +                  } else if (color[v] == ColorTraits::gray()) {
   1.471 +                    // find minimum flow on the cycle
   1.472 +                    FlowValue delta = residual_capacity[a];
   1.473 +                    while (1) {
   1.474 +                      BOOST_USING_STD_MIN();
   1.475 +                      delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]);
   1.476 +                      if (v == u)
   1.477 +                        break;
   1.478 +                      else
   1.479 +                        v = target(*current[v], g);
   1.480 +                    }
   1.481 +                    // remove delta flow units
   1.482 +                    v = u;
   1.483 +                    while (1) {
   1.484 +                      a = *current[v];
   1.485 +                      residual_capacity[a] -= delta;
   1.486 +                      residual_capacity[reverse_edge[a]] += delta;
   1.487 +                      v = target(a, g);
   1.488 +                      if (v == u)
   1.489 +                        break;
   1.490 +                    }
   1.491 +
   1.492 +                    // back-out of DFS to the first saturated edge
   1.493 +                    restart = u;
   1.494 +                    for (v = target(*current[u], g); v != u; v = target(a, g)){
   1.495 +                      a = *current[v];
   1.496 +                      if (color[v] == ColorTraits::white() 
   1.497 +                          || is_saturated(a)) {
   1.498 +                        color[target(*current[v], g)] = ColorTraits::white();
   1.499 +                        if (color[v] != ColorTraits::white())
   1.500 +                          restart = v;
   1.501 +                      }
   1.502 +                    }
   1.503 +                    if (restart != u) {
   1.504 +                      u = restart;
   1.505 +                      ++current[u];
   1.506 +                      break;
   1.507 +                    }
   1.508 +                  } // else if (color[v] == ColorTraits::gray())
   1.509 +                } // if (capacity[a] == 0 ...
   1.510 +              } // for out_edges(u, g)  (though "u" changes during loop)
   1.511 +              
   1.512 +              if (current[u] == out_edges(u, g).second) {
   1.513 +                // scan of i is complete
   1.514 +                color[u] = ColorTraits::black();
   1.515 +                if (u != src) {
   1.516 +                  if (bos_null) {
   1.517 +                    bos = u;
   1.518 +                    bos_null = false;
   1.519 +                    tos = u;
   1.520 +                  } else {
   1.521 +                    topo_next[u] = tos;
   1.522 +                    tos = u;
   1.523 +                  }
   1.524 +                }
   1.525 +                if (u != r) {
   1.526 +                  u = parent[u];
   1.527 +                  ++current[u];
   1.528 +                } else
   1.529 +                  break;
   1.530 +              }
   1.531 +            } // while (1)
   1.532 +          } // if (color[u] == white && excess_flow[u] > 0 & ...)
   1.533 +        } // for all vertices in g
   1.534 +
   1.535 +        // return excess flows
   1.536 +        // note that the sink is not on the stack
   1.537 +        if (! bos_null) {
   1.538 +          for (u = tos; u != bos; u = topo_next[u]) {
   1.539 +            ai = out_edges(u, g).first;
   1.540 +            while (excess_flow[u] > 0 && ai != out_edges(u, g).second) {
   1.541 +              if (capacity[*ai] == 0 && is_residual_edge(*ai))
   1.542 +                push_flow(*ai);
   1.543 +              ++ai;
   1.544 +            }
   1.545 +          }
   1.546 +          // do the bottom
   1.547 +          u = bos;
   1.548 +          ai = out_edges(u, g).first;
   1.549 +          while (excess_flow[u] > 0) {
   1.550 +            if (capacity[*ai] == 0 && is_residual_edge(*ai))
   1.551 +              push_flow(*ai);
   1.552 +            ++ai;
   1.553 +          }
   1.554 +        }
   1.555 +        
   1.556 +      } // convert_preflow_to_flow()
   1.557 +
   1.558 +      //=======================================================================
   1.559 +      inline bool is_flow()
   1.560 +      {
   1.561 +        vertex_iterator u_iter, u_end;
   1.562 +        out_edge_iterator ai, a_end;
   1.563 +
   1.564 +        // check edge flow values
   1.565 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.566 +          for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) {
   1.567 +            edge_descriptor a = *ai;
   1.568 +            if (capacity[a] > 0)
   1.569 +              if ((residual_capacity[a] + residual_capacity[reverse_edge[a]]
   1.570 +                   != capacity[a] + capacity[reverse_edge[a]])
   1.571 +                  || (residual_capacity[a] < 0)
   1.572 +                  || (residual_capacity[reverse_edge[a]] < 0))
   1.573 +              return false;
   1.574 +          }
   1.575 +        }
   1.576 +        
   1.577 +        // check conservation
   1.578 +        FlowValue sum;  
   1.579 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
   1.580 +          vertex_descriptor u = *u_iter;
   1.581 +          if (u != src && u != sink) {
   1.582 +            if (excess_flow[u] != 0)
   1.583 +              return false;
   1.584 +            sum = 0;
   1.585 +            for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) 
   1.586 +              if (capacity[*ai] > 0)
   1.587 +                sum -= capacity[*ai] - residual_capacity[*ai];
   1.588 +              else
   1.589 +                sum += residual_capacity[*ai];
   1.590 +
   1.591 +            if (excess_flow[u] != sum)
   1.592 +              return false;
   1.593 +          }
   1.594 +        }
   1.595 +
   1.596 +        return true;
   1.597 +      } // is_flow()
   1.598 +
   1.599 +      bool is_optimal() {
   1.600 +        // check if mincut is saturated...
   1.601 +        global_distance_update();
   1.602 +        return distance[src] >= n;
   1.603 +      }
   1.604 +
   1.605 +      void print_statistics(std::ostream& os) const {
   1.606 +        os << "pushes:     " << push_count << std::endl
   1.607 +           << "relabels:   " << relabel_count << std::endl
   1.608 +           << "updates:    " << update_count << std::endl
   1.609 +           << "gaps:       " << gap_count << std::endl
   1.610 +           << "gap nodes:  " << gap_node_count << std::endl
   1.611 +           << std::endl;
   1.612 +      }
   1.613 +
   1.614 +      void print_flow_values(std::ostream& os) const {
   1.615 +        os << "flow values" << std::endl;
   1.616 +        vertex_iterator u_iter, u_end;
   1.617 +        out_edge_iterator ei, e_end;
   1.618 +        for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
   1.619 +          for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)
   1.620 +            if (capacity[*ei] > 0)
   1.621 +              os << *u_iter << " " << target(*ei, g) << " " 
   1.622 +                 << (capacity[*ei] - residual_capacity[*ei]) << std::endl;
   1.623 +        os << std::endl;
   1.624 +      }
   1.625 +
   1.626 +      //=======================================================================
   1.627 +
   1.628 +      Graph& g;
   1.629 +      vertices_size_type n;
   1.630 +      vertices_size_type nm;
   1.631 +      EdgeCapacityMap capacity;
   1.632 +      vertex_descriptor src;
   1.633 +      vertex_descriptor sink;
   1.634 +      VertexIndexMap index;
   1.635 +
   1.636 +      // will need to use random_access_property_map with these
   1.637 +      std::vector< FlowValue > excess_flow;
   1.638 +      std::vector< out_edge_iterator > current;
   1.639 +      std::vector< distance_size_type > distance;
   1.640 +      std::vector< default_color_type > color;
   1.641 +
   1.642 +      // Edge Property Maps that must be interior to the graph
   1.643 +      ReverseEdgeMap reverse_edge;
   1.644 +      ResidualCapacityEdgeMap residual_capacity;
   1.645 +
   1.646 +      LayerArray layers;
   1.647 +      std::vector< list_iterator > layer_list_ptr;
   1.648 +      distance_size_type max_distance;  // maximal distance
   1.649 +      distance_size_type max_active;    // maximal distance with active node
   1.650 +      distance_size_type min_active;    // minimal distance with active node
   1.651 +      boost::queue<vertex_descriptor> Q;
   1.652 +
   1.653 +      // Statistics counters
   1.654 +      long push_count;
   1.655 +      long update_count;
   1.656 +      long relabel_count;
   1.657 +      long gap_count;
   1.658 +      long gap_node_count;
   1.659 +
   1.660 +      inline double global_update_frequency() { return 0.5; }
   1.661 +      inline vertices_size_type alpha() { return 6; }
   1.662 +      inline long beta() { return 12; }
   1.663 +
   1.664 +      long work_since_last_update;
   1.665 +    };
   1.666 +
   1.667 +  } // namespace detail
   1.668 +  
   1.669 +  template <class Graph, 
   1.670 +            class CapacityEdgeMap, class ResidualCapacityEdgeMap,
   1.671 +            class ReverseEdgeMap, class VertexIndexMap>
   1.672 +  typename property_traits<CapacityEdgeMap>::value_type
   1.673 +  push_relabel_max_flow
   1.674 +    (Graph& g, 
   1.675 +     typename graph_traits<Graph>::vertex_descriptor src,
   1.676 +     typename graph_traits<Graph>::vertex_descriptor sink,
   1.677 +     CapacityEdgeMap cap, ResidualCapacityEdgeMap res,
   1.678 +     ReverseEdgeMap rev, VertexIndexMap index_map)
   1.679 +  {
   1.680 +    typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;
   1.681 +    
   1.682 +    detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, 
   1.683 +      ReverseEdgeMap, VertexIndexMap, FlowValue>
   1.684 +      algo(g, cap, res, rev, src, sink, index_map);
   1.685 +    
   1.686 +    FlowValue flow = algo.maximum_preflow();
   1.687 +    
   1.688 +    algo.convert_preflow_to_flow();
   1.689 +    
   1.690 +    assert(algo.is_flow());
   1.691 +    assert(algo.is_optimal());
   1.692 +    
   1.693 +    return flow;
   1.694 +  } // push_relabel_max_flow()
   1.695 +  
   1.696 +  template <class Graph, class P, class T, class R>
   1.697 +  typename detail::edge_capacity_value<Graph, P, T, R>::type
   1.698 +  push_relabel_max_flow
   1.699 +    (Graph& g, 
   1.700 +     typename graph_traits<Graph>::vertex_descriptor src,
   1.701 +     typename graph_traits<Graph>::vertex_descriptor sink,
   1.702 +     const bgl_named_params<P, T, R>& params)
   1.703 +  {
   1.704 +    return push_relabel_max_flow
   1.705 +      (g, src, sink,
   1.706 +       choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
   1.707 +       choose_pmap(get_param(params, edge_residual_capacity), 
   1.708 +                   g, edge_residual_capacity),
   1.709 +       choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
   1.710 +       choose_const_pmap(get_param(params, vertex_index), g, vertex_index)
   1.711 +       );
   1.712 +  }
   1.713 +
   1.714 +  template <class Graph>
   1.715 +  typename property_traits<
   1.716 +    typename property_map<Graph, edge_capacity_t>::const_type
   1.717 +  >::value_type
   1.718 +  push_relabel_max_flow
   1.719 +    (Graph& g, 
   1.720 +     typename graph_traits<Graph>::vertex_descriptor src,
   1.721 +     typename graph_traits<Graph>::vertex_descriptor sink)
   1.722 +  {
   1.723 +    bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
   1.724 +    return push_relabel_max_flow(g, src, sink, params);
   1.725 +  }
   1.726 +
   1.727 +} // namespace boost
   1.728 +
   1.729 +#endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP
   1.730 +