From 8b11baa5cf1fa602952c97e353ac41fa01e36356 Mon Sep 17 00:00:00 2001 From: Jeremy Siek Date: Sun, 12 Nov 2000 23:42:35 +0000 Subject: [PATCH] making some changes with regards to the hi_pr.c implementation [SVN r8198] --- include/boost/graph/maximum_flow.hpp | 371 ++++++++++++++------------- 1 file changed, 197 insertions(+), 174 deletions(-) diff --git a/include/boost/graph/maximum_flow.hpp b/include/boost/graph/maximum_flow.hpp index a80dadda..ad1d14bd 100644 --- a/include/boost/graph/maximum_flow.hpp +++ b/include/boost/graph/maximum_flow.hpp @@ -7,7 +7,7 @@ namespace boost { - enum edge_sister_t { edge_sister }; + enum edge_reverse_t { edge_reverse }; enum edge_residual_capacity_t { edge_residual_capacity }; namespace detail { @@ -17,7 +17,7 @@ namespace boost { "On Implementing Push-Relabel Method for the Maximum Flow Problem" by B.V. Cherkassky and A.V. Goldberg, IPCO '95, 157--171. - Also based on the h_prf.c code written by the above authors. + Also based on the h_prf.c and hi_pr.c code written by the above authors. This implements the highest-label version of the push-relabel method with the global relabeling and gap relabeling heuristics. @@ -58,25 +58,26 @@ namespace boost { template struct preflow_layer { preflow_layer() - : push_list(new list_node()), - transit_list(new list_node()) + : active_vertices(new list_node()), + inactive_vertices(new list_node()) { } ~preflow_layer() { - dlist_clear(push_list, delete_node()); - delete push_list; - dlist_clear(transit_list, delete_node()); - delete transit_list; + dlist_clear(active_vertices, delete_node()); + delete active_vertices; + dlist_clear(inactive_vertices, delete_node()); + delete inactive_vertices; } typedef list_node list_node; - list_node* push_list; // nodes with positive excess - list_node* transit_list; // nodes with zero excess - - // push_list and transit_list point to the dummy node of the dlist. - // next(push_list) is the first node of the push list. - // next(transit_list) is the first node of the transit list. + list_node* active_vertices; // nodes with positive excess + list_node* inactive_vertices; // nodes with zero excess + + // active_vertices and inactive_vertices point to the dummy node + // of the dlist. next(active_vertices) is the first node of the + // push list. next(inactive_vertices) is the first node of the + // transit list. }; - // Graph must have edge_sister and edge_rcap properties + // Graph must have edge_reverse and edge_rcap properties template integer @@ -91,19 +92,45 @@ namespace boost { typedef typename Traits::out_edge_iterator out_edge_iterator; typedef typename Traits::vertices_size_type vertices_size_type; - typedef typename property_map::type SisterMap; - typedef typename property_map::type + typedef typename property_map::type SisterMap; + typedef typename property_map::type ResidualCapacityMap; typedef preflow_layer Layer; typedef std::vector< Layer > LayerArray; typedef typename LayerArray::iterator layer_iterator; - typedef typename LayerArray::size_type rank_size_type; + typedef typename LayerArray::size_type distance_size_type; typedef color_traits ColorTraits; typedef list_node list_node; + + //======================================================================= + // Layer List Management Functions + + void add_active(vertex_descriptor u, Layer& layer) { + list_node* u_node = new list_node(u); + dlist_insert_before(next(layer.active_vertices), u_node, next, prev); + max_active = std::max(distance[j], max_active); + min_active = std::min(distance[j], min_active); + layer_list_ptr[u] = u_node; + } + void remove_active(vertex_descriptor u) { + dlist_remove(layer_list_ptr[u], next, prev); + delete layer_list_ptr[u]; + } + + void add_inactive(vertex_descriptor u, Layer& layer) { + list_node* u_node = new list_node(u); + dlist_insert_before(next(layer.inactive_vertices), u_node, next, prev); + layer_list_ptr[u] = u_node; + } + void remove_inactive(vertex_descriptor u) { + dlist_remove(layer_list_ptr[u], next, prev); + delete layer_list_ptr[u]; + } + //======================================================================= push_relabel(Graph& g_, EdgeCapacityMap cap_, @@ -115,9 +142,9 @@ namespace boost { excess_flow(num_vertices(g_)), layer_list_ptr(num_vertices(g_)), current(num_vertices(g_)), - rank(num_vertices(g_)), + distance(num_vertices(g_)), color(num_vertices(g_)), - sister(get(edge_sister, g_)), + reverse_edge(get(edge_sister, g_)), residual_capacity(get(edge_residual_capacity, g_)) { vertex_iterator i_iter, end; @@ -125,7 +152,7 @@ namespace boost { excess_flow[*i_iter] = 0; excess_flow[src] = std::numeric_limits::max(); - lmax = num_vertices(g) - 1; + max_distance = num_vertices(g) - 1; } //======================================================================= @@ -133,95 +160,101 @@ namespace boost { // (well, actually the reverse of the residual graph). // Would be cool to have a graph view adaptor for hiding certain // edges, like the saturated (non-residual) edges in this case. - // Goldberg's implementation abused "rank" for the coloring. - void global_rank_update() + // Goldberg's implementation abused "distance" for the coloring. + void global_distance_update() { vertex_iterator i_iter, i_end; for (tie(i_iter,i_end) = vertices(g); i_iter != i_end; ++i_iter) { color[*i_iter] = ColorTraits::white(); + distance[*i_iter] = n; + } color[sink] = ColorTraits::gray(); + distance[sink] = 0; + + for (distance_size_type l = 0; l != max_distance; ++l) { + dlist_clear(layers[l].active_vertices); + dlist_clear(layers[l].inactive_vertices); + } + + max_distance = max_active = 0; + min_active = n; Q.push(sink); - lmax = lmax_push = 0; - lmin_push = n; - while (! Q.empty()) { vertex_descriptor i = Q.top(); Q.pop(); - rank_size_type j_rank = rank[i] + 1; + distance_size_type j_distance = distance[i] + 1; out_edge_iterator ai, a_end; for (tie(ai, a_end) = out_edges(i, g); ai != a_end; ++ai) { edge_descriptor a = *ai; vertex_descriptor j = target(a, g); if ( color[j] == ColorTraits::white() - && residual_capacity[sister[a]] > 0 ) { - rank[j] = j_rank; + && residual_capacity[reverse_edge[a]] > 0 ) { + distance[j] = j_distance; color[j] = ColorTraits::gray(); current[j] = out_edges(j, g).first; - Layer& layer = layers[j_rank]; - lmax = std::max(j_rank, lmax); + max_distance = std::max(j_distance, max_distance); + + if (excess_flow[j] > 0) + add_active(j, layers[j_distance]); + else + add_inactive(j, layers[j_distance]); - if (excess_flow[j] > 0) { - list_node* pj = new list_node(j); - dlist_insert_before(next(layer.push_list), pj, next, prev); - layer_list_ptr[j] = pj; - lmax_push = std::max(j_rank, lmax_push); - lmin_push = std::min(j_rank, lmin_push); - } else { - list_node* pj = new list_node(j); - dlist_insert_before(next(layer.transit_list), pj, next, prev); - layer_list_ptr[j] = pj; - } Q.push(j); - } // if (rank[j] == n && ...) - } // for out_edges(i) - } // while (! Q.empty()) - } // global_rank_update() + } + } + } + } // global_distance_update() //======================================================================= - // This function is called "push" in Goldberg's implementation, - // but it is called "discharge" in the paper. - bool discharge(vertex_descriptor i) + // This function is called "push" in Goldberg's h_prf implementation, + // but it is called "discharge" in the paper and in hi_pr.c. + void discharge(vertex_descriptor i) { - rank_size_type next_layer_rank = rank[i] - 1; + assert(excess_flow[i] > 0); + while (1) { + distance_size_type next_layer_distance = distance[i] - 1; + out_edge_iterator ai, ai_end; + for (ai = current[i], ai_end = out_edges(i, g).second; + ai != ai_end; ++ai) { + edge_descriptor a = *ai; + if (residual_capacity[a] > 0) { + vertex_descriptor j = target(a, g); + if (distance[j] == next_layer_distance) { + if (j != sink && excess_flow[j] == 0) { + remove_inactive(j); + add_active(j, layers[next_layer_distance]); + } + push_flow(a); + if (excess_flow[i] == 0) + break; + } + } + } // for out_edges of i starting from current - out_edge_iterator ai, a_end; - ai = current[i]; a_end = out_edges(i, g).second; - for (; ai != a_end; ++ai) { - edge_descriptor a = *ai; - vertex_descriptor j = target(a, g); + Layer& layer = layers[distance[i]]; - if (rank[j] == next_layer_rank) { // if j belongs to the next layer - - if (next_layer_rank > 0) { - next_layer = layers[next_layer_rank]; - - if (excess_flow[j] == 0) { // j already had zero excess flow - dlist_remove(layer_list_ptr[j], next, prev); - dlist_insert_before(next(next_layer.push_list), - layer_list_ptr[j], next, prev); - lmin_push = std::min(next_layer_rank, lmin_push); - } // if (excess_flow[j] == 0) - - } // if (next_layer_rank > 0) - - push(a); - - if (excess_flow[i] == 0) + if (ai == ai_end) { // i must be relabeled + relabel(i); + if (dlist_empty(layer.active_vertices) + && dlist_empty(layer.inactive_vertices)) + gap(distance[i]); + if (distance[i] == n) break; + } else { // i is no longer active + current[i] = ai; + add_inactive(i, layer); + break; + } - } // if (i,j) is admissible - } // for out_edges of i from current - - current[i] = ai; - return ai == out_edges(i,g).second ? true : false; // end of list? + } // while (1) } // discharge() //======================================================================= // This corresponds to the "push" update operation of the paper, // not the "push" function in Goldberg's implementation. - void push(edge_descriptor u_v) + void push_flow(edge_descriptor u_v) { vertex_descriptor u = source(u_v, g), @@ -231,74 +264,58 @@ namespace boost { = std::min(excess_flow[u], residual_capacity[u_v]); residual_capacity[u_v] -= flow_delta; - residual_capacity[sister[u_v]] += flow_delta; + residual_capacity[reverse_edge[u_v]] += flow_delta; excess_flow[u] -= flow_delta; excess_flow[v] += flow_delta; - } + } // push_flow() //======================================================================= - rank_size_type relabel(vertex_descriptor i) + distance_size_type relabel(vertex_descriptor i) { - rank_size_type j_rank = num_vertices(g); - rank[i] = j_rank; + distance_size_type min_distance = num_vertices(g); + distance[i] = min_distance; // Examine the residual out-edges of vertex i, choosing the - // edge whose target vertex has the minimal rank. + // edge whose target vertex has the minimal distance. // This is the highest-label or "HL" optimization. - out_edge_iterator ai, a_end, a_j; + out_edge_iterator ai, a_end, min_edge_iter; for (tie(ai, a_end) = out_edges(i, g); ai != a_end; ++ai) { edge_descriptor a = *ai; vertex_descriptor j = target(a, g); - if (residual_capacity[a] > 0 && rank[j] < j_rank) { - j_rank = rank[j]; - a_j = ai; + if (residual_capacity[a] > 0 && distance[j] < min_distance) { + min_distance = distance[j]; + min_edge_iter = ai; } } // for all out edges of i - ++j_rank; - if (j_rank < n) { - rank[i] = j_rank; // this is the main action - current[i] = a_j; - Layer& layer = layers[j_rank]; - if (excess_flow[i] > 0) { - list_node* pi = new list_node(i); - dlist_insert_before(next(layer.push_list), pi, next, prev); - layer_list_ptr[i] = pi; - lmax_push = std::max(j_rank, lmax_push); - lmin_mush = std::min(j_rank, lmin_push); - } else { - list_node* pi = new list_node(i); - dlist_insert_before(next(layer.transit_list), pi, next, prev); - layer_list_ptr[i] = pi; - } // if (excess_flow[i] > 0) - lmax = std::max(j_rank, lmax); - } // (j_rank < n) + ++min_distance; + if (min_distance < n) { + distance[i] = min_distance; // this is the main action + current[i] = min_edge_iter; + max_distance = std::max(min_distance, max_distance); + } // (min_distance < n) - return j_rank; + return min_distance; } // relabel() //======================================================================= // cleanup beyond the gap - void gap(layer_iterator empty_layer) + void gap(distance_size_type empty_distance) { - rank_size_type r; // rank of layer before the current layer - r = (empty_layer - layers.begin()) - 1; + distance_size_type r; // distance of layer before the current layer + r = empty_distance - 1; - for (layer_iterator l = empty_layer + 1; - l != layers.begin() + lmax; ++l) { + for (layer_iterator l = layers.begin() + empty_distance + 1; + l != layers.begin() + max_distance; ++l) { list_node* i; - for (i = next(l->push_list); i != l->push_list; i = next(i)) - rank[i->m_data] = n; - - for (i = next(l->transit_list); i != l->transit_list; i = next(i)) - rank[i->m_data] = n; - - dlist_clear(l->push_list); - dlist_clear(l->trans_list); + for (i = next(l->inactive_vertices); + i != l->inactive_vertices; i = next(i)) + distance[i->m_data] = n; + dlist_clear(l->inactive_vertices); } - lmax = r; - lmax_push = r; + max_distance = r; + max_active = r; } //======================================================================= @@ -309,7 +326,7 @@ namespace boost { // Return the nodes with excess flow in topological order. // // Unlike the prefl_to_flow() implementation, we use - // "color" instead of "rank" for the DFS labels + // "color" instead of "distance" for the DFS labels // "parent" instead of nl_prev for the DFS tree // "topo_next" instead of nl_next for the topological ordering void convert_preflow_to_flow() @@ -317,8 +334,10 @@ namespace boost { vertex_iterator i_iter, i_end; out_edge_iterator ai, a_end; - vertex_descriptor r, restart; - vertex_iterator tos, bos; + vertex_descriptor r, restart, i; + + vertex_descriptor tos, bos; + bool bos_null = true; std::vector parent(n); std::vector topo_next(n); @@ -336,7 +355,7 @@ namespace boost { parent[i] = i; current[i] = out_edges(i, g).first; } - + // eliminate flow cycles and topologically order the vertices for (tie(i_iter, i_end) = vertices(g); i_iter != i_end; ++i_iter) { i = *i_iter; if (color[i] == ColorTraits::white() @@ -347,8 +366,7 @@ namespace boost { while (1) { for (; current[i] != out_edges(i, g).second; ++current[i]) { a = current[i]; - if ( capacity[a] == 0 && residual_capacity[a] > 0 - && target(a, g) != src && target(a, g) != sink ) { + if ( capacity[a] == 0 && residual_capacity[a]) { j = target(a, g); if (color[j] == ColorTraits::white()) { color[j] = ColorTraits::gray(); @@ -370,13 +388,13 @@ namespace boost { while (1) { a = *current[j]; residual_capacity[a] -= delta; - residual_capacity[sister[a]] += delta; + residual_capacity[reverse_edge[a]] += delta; j = target(a, g); if (j == i) break; } - // back-out DFS to the first zeroed edge + // back-out of DFS to the first saturated edge restart = i; for (j = target(*current[i], g); j != i; j = target(a, g)){ a = current[j]; @@ -395,49 +413,54 @@ namespace boost { } // else if (color[j] == ColorTraits::gray()) } // if (capacity[a] == 0 ... } // for out_edges(i, g) (though "i" changes during loop) - + if (current[i] == out_edges(i, g).second) { // scan of i is complete color[i] = ColorTraits::black(); if (i != src) { - if (bos == NULL) { - bos = i_iter; - tos = i_iter; + if (bos_null) { + bos = i; + bos_null = false; + tos = i; } else { - topo_next[i] = *tos; - tos = i_iter; + topo_next[i] = tos; + tos = i; } } if (i != r) { - i = topo_next[i]; - ++current[i]; + i = topo_next[i]; + ++current[i]; } else break; } } // while (1) - } // if color == white ... - - i_iter = vertices(g).first + i; + // In case i was changed, adjust i_iter to match + i_iter = vertices(g).first + get(index, i); + } // if (color[i] == white && excess_flow[i] > 0 & ...) } // for all vertices in g - if (bos != vertices(g).second) { - i_iter = tos; - do { - i = *i_iter; + // return excesses + // note that sink is not on the stack + if (! bos_null) { + for (i = tos; i != bos; i = topo_next[i]) { ai = out_edges(i, g).first; a = *ai; - while (excess_flow[i] > 0 && ai != out_edges(i, g).second) { + while (excess_flow[i] > 0) { if (capacity[a] == 0 && residual_capacity[a] > 0) - push(a); + push_flow(a); ++ai; } - if (i == *bos) - break; - else - i = topo_next[i]; - } while ( 1 ); + } + // now at the bottom + i = bos; + a = out_edge(i, g).first; + while (excess_flow[i] > 0) { + if (capacity[a] == 0 && residual_capacity[a] > 0) + push_flow(a); + ++ai; + } } - + } // convert_preflow_to_flow() //======================================================================= @@ -453,20 +476,20 @@ namespace boost { std::vector< FlowValue > excess_flow; std::vector< list_node* > layer_list_ptr; std::vector< out_edge_iterator > current; - std::vector< rank_size_type > rank; + std::vector< distance_size_type > distance; std::vector color; // Edge Property Maps that must be interior to the graph - SisterMap sister; + SisterMap reverse_edge; ResidualCapacityMap residual_capacity; LayerArray layers; - rank_size_type lmax; // maximal layer - rank_size_type lmax_push; // maximal layer with excess node - rank_size_type lmin_push; // minimal layer with excess node + distance_size_type max_distance; // maximal distance + distance_size_type max_active; // maximal distance with active node + distance_size_type min_active; // minimal distance with active node std::queue Q; - inline double global_update_frequency() { return 1.0; } + inline double global_update_frequency() { return 0.5; } list_node_next next; list_node_prev prev; @@ -486,30 +509,30 @@ namespace boost { detail::push_relabel algo(g, cap, src, sink, index_map); - algo.global_rank_update(); + algo.global_distance_update(); typename graph_traits::vertices_size_type num_relabels = 0, n = num_vertices(g); - while (algo.lmax_push >= algo.lmin_push) { + while (algo.max_active >= algo.min_active) { - Layer& layer = algo.layers[algo.lmax_push]; + Layer& layer = algo.layers[algo.max_active]; - list_node* i_node = next(layer.push_list); - if (i_node == layer.push_list) - --algo.lmax_push; + list_node* i_node = next(layer.active_vertices); + if (i_node == layer.active_vertices) + --algo.max_active; else { - dlist_remove(next(layer.push_list), next, prev); + dlist_remove(next(layer.active_vertices), next, prev); if (algo.discharge(*i_node)) { // returns true if done with out-edges algo.relabel(*i_node); ++num_relabels; // Gap relabelling heuristic - if (dlist_empty(layer.push_list) && dlist_empty(layer.transit_list)) + if (dlist_empty(layer.active_vertices) && dlist_empty(layer.inactive_vertices)) algo.gap(layer); // Global relabelling heuristic if (num_relabels > algo.global_update_frequency() * n) { - algo.global_rank_update(); + algo.global_distance_update(); num_relabels = 0; } @@ -519,7 +542,7 @@ namespace boost { } } - } // while (algo.lmax_push >= algo.lmin_push) + } // while (algo.max_active >= algo.min_active) flow += algo.excess_flow[sink]; algo.convert_preflow_to_flow(); @@ -588,20 +611,20 @@ namespace boost { // Vertex Properties: // - // layer_list_ptr[v] points into push_list or trans_list + // layer_list_ptr[v] points into active_vertices or inactive_vertices // current[v] (value_type == out_edge_iterator) - // rank[v] also known as label + // distance[v] also known as label // excess_flow[v] // Edge Properties: // - // sister[(u,v)] => (v,u) + // reverse_edge[(u,v)] => (v,u) // residual_capacity[(u,v)] #if 0 -typedef std::list RankSet; -typedef std::vector ArrayOfRankSets; +typedef std::list DistanceSet; +typedef std::vector ArrayOfDistanceSets; ArrayOfLabelSets B; ArrayOfLabelSets::size_type b;