From d938af040478198c67659ebf95052fd3c3d2eb86 Mon Sep 17 00:00:00 2001 From: Jeremy Siek Date: Sun, 12 Nov 2000 20:28:37 +0000 Subject: [PATCH] more work on max-flow, getting closer to finishing [SVN r8179] --- include/boost/graph/detail/list_base.hpp | 209 ++++++++++++++ include/boost/graph/maximum_flow.hpp | 341 ++++++++++++++--------- 2 files changed, 426 insertions(+), 124 deletions(-) create mode 100644 include/boost/graph/detail/list_base.hpp diff --git a/include/boost/graph/detail/list_base.hpp b/include/boost/graph/detail/list_base.hpp new file mode 100644 index 00000000..6cb166c0 --- /dev/null +++ b/include/boost/graph/detail/list_base.hpp @@ -0,0 +1,209 @@ +#ifndef BOOST_LIST_BASE_HPP +#define BOOST_LIST_BASE_HPP + +// Perhaps this should go through formal review, and move to . + +/* + An alternate interface idea: + Extend the std::list functionality by creating remove/insert + functions that do not require the container object! + */ + +namespace boost { + namespace detail { + + //========================================================================= + // Linked-List Generic Implementation Functions + + template + inline Node + slist_insert_after(Node pos, Node x, + Next next) + { + next(x) = next(pos); + next(pos) = x; + return x; + } + + // return next(pos) or next(next(pos)) ? + template + inline Node + slist_remove_after(Node pos, + Next next) + { + Node n = next(pos); + next(pos) = next(n); + return n; + } + + template + inline Node + slist_remove_range(Node before_first, Node last, + Next next) + { + next(before_first) = last; + return last; + } + + template + inline Node + slist_previous(Node head, Node x, Node nil, + Next next) + { + while (head != nil && next(head) != x) + head = next(head); + return head; + } + + template + inline void + slist_splice_after(Node pos, Node before_first, Node before_last, + Next next) + { + if (pos != before_first && pos != before_last) { + Node first = next(before_first); + Node after = next(pos); + next(before_first) = next(before_last); + next(pos) = first; + next(before_last) = after; + } + } + + template + inline Node + slist_reverse(Node node, Node nil, + Next next) + { + Node result = node; + node = next(node); + next(result) = nil; + while(node) { + Node next = next(node); + next(node) = result; + result = node; + node = next; + } + return result; + } + + template + inline std::size_t + slist_size(Node head, Node nil, + Next next) + { + std::size_t s = 0; + for ( ; head != nil; head = next(head)) + ++s; + return s; + } + + template + class slist_iterator_policies + { + public: + explicit slist_iterator_policies(const Next& n, const Data& d) + : m_next(n), m_data(d) { } + + template + Reference dereference(type, const Node& x) const + { return m_data(x); } + + template + void increment(Node& x) const + { x = m_next(x); } + + template + bool equal(Node& x, Node& y) const + { return x == y; } + + protected: + Next m_next; + Data m_data; + }; + + //=========================================================================== + // Doubly-Linked List Generic Implementation Functions + + template + inline void + dlist_insert_before(Node pos, Node x, + Next next, Prev prev) + { + next(x) = pos; + prev(x) = prev(pos); + next(prev(pos)) = x; + prev(pos) = x; + } + + template + void + dlist_remove(Node pos, + Next next, Prev prev) + { + Node next_node = next(pos); + Node prev_node = prev(pos); + next(prev_node) = next_node; + prev(next_node) = prev_node; + } + + // This deletes every node in the list except the + // "dummy" node. + template + inline void + dlist_clear(Node dummy, Delete del) + { + Node i, tmp; + i = next(head); + while (i != head) { + tmp = i; + i = next(i) + del(tmp); + } + } + + template + inline bool + dlist_empty(Node dummy) + { + return next(dummy) == dummy; + } + + template + void + dlist_transfer(Node pos, Node first, Node last, + Next next, Prev prev) + { + if (pos != last) { + // Remove [first,last) from its old position + next(prev(last)) = pos; + next(prev(first)) = last; + next(prev(pos)) = first; + + // Splice [first,last) into its new position + Node tmp = prev(pos); + prev(pos) = prev(last); + prev(last) = prev(first); + prev(first) = tmp; + } + } + + template + class dlist_iterator_policies + : public slist_iterator_policies + { + typedef slist_iterator_policies Base; + public: + template + void decrement(Node& x) const + { x = m_prev(x); } + + dlist_iterator_policies(Next n, Prev p, Data d) + : Base(n,d), m_prev(p) { } + protected: + Prev m_prev; + }; + + } // namespace detail +} // namespace boost + +#endif // BOOST_LIST_BASE_HPP diff --git a/include/boost/graph/maximum_flow.hpp b/include/boost/graph/maximum_flow.hpp index f972791a..a80dadda 100644 --- a/include/boost/graph/maximum_flow.hpp +++ b/include/boost/graph/maximum_flow.hpp @@ -3,10 +3,12 @@ // UNDER CONSTRUCTION +#include + namespace boost { enum edge_sister_t { edge_sister }; - enum edge_rcap_t { edge_sister }; + enum edge_residual_capacity_t { edge_residual_capacity }; namespace detail { @@ -15,26 +17,72 @@ 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. + This implements the highest-label version of the push-relabel method with the global relabeling and gap relabeling heuristics. need to write a DIMACS graph format reader... + + The "rank" or layers used here and in Goldberg's implementation + is the reverse of the distance discussed in the paper. + + Each layer consists of all active nodes with the same rank. + An active node has positive excess flow and its + distance is less than n (it is not blocked). */ + template + struct list_node { + explicit list_node(const T& x) + : m_data(x), m_next(*this), m_prev(*this) { } + T m_data; + list_node* m_next; + list_node* m_prev; + }; + template + struct list_node_next { + T& operator()(list_node* x) { return x->m_next; } + const T& operator()(const list_node* x) { return x->m_next; } + }; + template + struct list_node_prev { + T& operator()(list_node* x) { return x->m_prev; } + const T& operator()(const list_node* x) { return x->m_prev; } + }; + template + struct delete_node { + void operator()(const list_node* x) { delete x; } + }; + template struct preflow_layer { - typedef std::list list_type; - typedef typename layer_list_type::iterator list_iterator; - LayerList push_list; // nodes with positive excess - LayerList trans_list; // nodes with zero excess + preflow_layer() + : push_list(new list_node()), + transit_list(new list_node()) + { } + ~preflow_layer() { + dlist_clear(push_list, delete_node()); + delete push_list; + dlist_clear(transit_list, delete_node()); + delete transit_list; + } + 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. }; // Graph must have edge_sister and edge_rcap properties template integer class FlowValue> - class preflow_push { + class push_relabel + { public: typedef graph_traits Traits; typedef typename Traits::vertex_descriptor vertex_descriptor; @@ -52,32 +100,46 @@ namespace boost { typedef typename LayerArray::iterator layer_iterator; typedef typename LayerArray::size_type rank_size_type; + typedef color_traits ColorTraits; + + typedef list_node list_node; + //======================================================================= - preflow_push(Graph& g_, + push_relabel(Graph& g_, EdgeCapacityMap cap_, vertex_descriptor src_, vertex_descriptor sink_, VertexIndexMap idx_) : g(g_), n(num_vertices(g_)), capacity(cap_), src(src_), sink(sink_), index(idx_), - sister(get(edge_sister, g_)), - residual_capacity(get(edge_rcap, g_)) + excess_flow(num_vertices(g_)), + layer_list_ptr(num_vertices(g_)), + current(num_vertices(g_)), + rank(num_vertices(g_)), + color(num_vertices(g_)), + sister(get(edge_sister, g_)), + residual_capacity(get(edge_residual_capacity, g_)) { vertex_iterator i_iter, end; for (tie(i_iter, end) = vertices(g); i_iter != end; ++i_iter) - excess_flow[*i_iter] = FlowValue(); + excess_flow[*i_iter] = 0; excess_flow[src] = std::numeric_limits::max(); lmax = num_vertices(g) - 1; } //======================================================================= + // This is a breadth-first search over the residual graph + // (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() { vertex_iterator i_iter, i_end; for (tie(i_iter,i_end) = vertices(g); i_iter != i_end; ++i_iter) { - rank[*i_iter] = n; - rank[sink] = 0; + color[*i_iter] = ColorTraits::white(); + color[sink] = ColorTraits::gray(); Q.push(sink); lmax = lmax_push = 0; @@ -92,20 +154,24 @@ namespace boost { for (tie(ai, a_end) = out_edges(i, g); ai != a_end; ++ai) { edge_descriptor a = *ai; vertex_descriptor j = target(a, g); - if ( rank[j] == n && residual_capacity[sister[a]] > 0 ) { + if ( color[j] == ColorTraits::white() + && residual_capacity[sister[a]] > 0 ) { rank[j] = j_rank; + color[j] = ColorTraits::gray(); current[j] = out_edges(j, g).first; Layer& layer = layers[j_rank]; lmax = std::max(j_rank, lmax); if (excess_flow[j] > 0) { - layer.push_list.push_front(j); - layer_list_iter[j] = layer.push_list.begin(); + 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 { - layer.transit_list.push_front(j); - layer_list_iter[j] = layer.transit_list.begin(); + 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 && ...) @@ -114,7 +180,9 @@ namespace boost { } // global_rank_update() //======================================================================= - bool push(vertex_descriptor i) + // This function is called "push" in Goldberg's implementation, + // but it is called "discharge" in the paper. + bool discharge(vertex_descriptor i) { rank_size_type next_layer_rank = rank[i] - 1; @@ -124,34 +192,50 @@ namespace boost { edge_descriptor a = *ai; vertex_descriptor j = target(a, g); - if (rank[j] == next_layer_rank) { + if (rank[j] == next_layer_rank) { // if j belongs to the next layer - flow_delta = std::min(excess_flow[i], residual_capacity[a]); - - residual_capacity[a] -= flow_delta; - residual_capacity[sister[a]] += flow_delta; - - if (next_layer_rank > 0) { + if (next_layer_rank > 0) { next_layer = layers[next_layer_rank]; - if (excess_flow[j] == 0) { - next_layer.transit_list.erase( layer_list_iter[j] ); - next_layer.push_list.push_front(j); + + 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) - excess_flow[j] += flow_delta; - excess_flow[i] -= flow_delta; + push(a); - if (excess_flow[i] == 0) + if (excess_flow[i] == 0) break; - } // if (rank[j] == next_layer_rank) + } // if (i,j) is admissible } // for out_edges of i from current current[i] = ai; - return ai == out_edges(i,g).second ? true : false; - } // push() + return ai == out_edges(i,g).second ? true : false; // end of list? + } // 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) + { + vertex_descriptor + u = source(u_v, g), + v = target(u_v, g); + + FlowValue flow_delta + = std::min(excess_flow[u], residual_capacity[u_v]); + + residual_capacity[u_v] -= flow_delta; + residual_capacity[sister[u_v]] += flow_delta; + + excess_flow[u] -= flow_delta; + excess_flow[v] += flow_delta; + } //======================================================================= rank_size_type relabel(vertex_descriptor i) @@ -159,31 +243,34 @@ namespace boost { rank_size_type j_rank = num_vertices(g); rank[i] = j_rank; + // Examine the residual out-edges of vertex i, choosing the + // edge whose target vertex has the minimal rank. + // This is the highest-label or "HL" optimization. out_edge_iterator ai, a_end, a_j; for (tie(ai, a_end) = out_edges(i, g); ai != a_end; ++ai) { edge_descriptor a = *ai; - if (residual_capacity[a] > 0) { - vertex_descriptor j = target(a, g); - if (rank[j] < j_rank) { + vertex_descriptor j = target(a, g); + if (residual_capacity[a] > 0 && rank[j] < j_rank) { j_rank = rank[j]; a_j = ai; - } } } // for all out edges of i ++j_rank; if (j_rank < n) { - rank[i] = j_rank; + rank[i] = j_rank; // this is the main action current[i] = a_j; Layer& layer = layers[j_rank]; if (excess_flow[i] > 0) { - layer.push_list.push_front(, i); - layer_list_iter[i] = layer.push_list.begin(); + 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 { - layer.transit_list.push_front(i); - layer_list_iter[i] = layer.transit_list.begin(); + 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) @@ -200,16 +287,15 @@ namespace boost { for (layer_iterator l = empty_layer + 1; l != layers.begin() + lmax; ++l) { - for (typename Layer::list_iterator i = l->push_list.begin(); - i != l->push_list.end(); ++i) - rank[*i] = n; + list_node* i; + for (i = next(l->push_list); i != l->push_list; i = next(i)) + rank[i->m_data] = n; - for (typename Layer::list_iterator i = l->trans_list.begin(); - i != l->push_list.end(); ++i) - rank[*i] = n; + for (i = next(l->transit_list); i != l->transit_list; i = next(i)) + rank[i->m_data] = n; - l->push_list.clear(); - l->trans_list.clear(); + dlist_clear(l->push_list); + dlist_clear(l->trans_list); } lmax = r; lmax_push = r; @@ -225,19 +311,17 @@ namespace boost { // Unlike the prefl_to_flow() implementation, we use // "color" instead of "rank" for the DFS labels // "parent" instead of nl_prev for the DFS tree - // "topo_order" instead of nl_next for the topological ordering + // "topo_next" instead of nl_next for the topological ordering void convert_preflow_to_flow() { vertex_iterator i_iter, i_end; out_edge_iterator ai, a_end; - // tos, bos, restart, r: vertex_iterator or vertex_descriptor? + vertex_descriptor r, restart; + vertex_iterator tos, bos; - std::vector color(n); - typedef color_traits ColorTr; - - std::vector parent(n); - std::vector topo_order; + std::vector parent(n); + std::vector topo_next(n); // handle self-loops for (tie(i_iter, i_end) = vertices(g); i_iter != i_end; ++i_iter) @@ -248,30 +332,31 @@ namespace boost { // initialize for (tie(i_iter, i_end) = vertices(g); i_iter != i_end; ++i_iter) { i = *i_iter; - color[i] = ColorTr::white(); + color[i] = ColorTraits::white(); + parent[i] = i; current[i] = out_edges(i, g).first; - parent[i] = i; } - for (tie(i_iter, i_end) = vertices(g); i_iter != i_end; ++i) { + for (tie(i_iter, i_end) = vertices(g); i_iter != i_end; ++i_iter) { i = *i_iter; - if ( color[i] == ColorTr::white() && excess_flow[i] > 0 - && i != src && i != sink ) { + if (color[i] == ColorTraits::white() + && excess_flow[i] > 0 + && i != src && i != sink ) { r = i; - color[r] = ColorTr::gray(); - while (1) { + color[r] = ColorTraits::gray(); + 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 ) { j = target(a, g); - if (color[j] == ColorTr::white()) { - color[j] = ColorTr::gray(); - parent[j] = i; - i = j; // ?? + if (color[j] == ColorTraits::white()) { + color[j] = ColorTraits::gray(); + parent[j] = i; + i = j; break; - } else if (color[j] == ColorTr::gray()) { - // find minimum flow on the cycle + } else if (color[j] == ColorTraits::gray()) { + // find minimum flow on the cycle delta = residual_capacity[a]; while (1) { delta = std::min(delta, residual_capacity[*current[j]]); @@ -280,25 +365,25 @@ namespace boost { else j = target(*current[j], g); } - // remove delta flow units - j = i; - while (1) { - a = *current[j]; - residual_capacity[a] -= delta; - residual_capacity[sister[a]] += delta; - j = target(a, g); - if (j == i) - break; - } + // remove delta flow units + j = i; + while (1) { + a = *current[j]; + residual_capacity[a] -= delta; + residual_capacity[sister[a]] += delta; + j = target(a, g); + if (j == i) + break; + } // back-out DFS to the first zeroed edge restart = i; - for (j = target(*current[i], g); j != i; j = target(a, g)) { + for (j = target(*current[i], g); j != i; j = target(a, g)){ a = current[j]; - if (color[j] == ColorTr::white() || - residual_capacity[a] == 0) { - color[target(*current[j], g)] = ColorTr::white(); - if (color[j] != ColorTr::white()) + if (color[j] == ColorTraits::white() || + residual_capacity[a] == 0) { + color[target(*current[j], g)] = ColorTraits::white(); + if (color[j] != ColorTraits::white()) restart = j; } } // for @@ -307,54 +392,52 @@ namespace boost { ++current[i]; break; } - } // else if (color[j] == ColorTr::gray()) + } // else if (color[j] == ColorTraits::gray()) } // if (capacity[a] == 0 ... - } // for + } // 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] = ColorTr::black(); + // scan of i is complete + color[i] = ColorTraits::black(); if (i != src) { if (bos == NULL) { bos = i_iter; tos = i_iter; } else { - layer_list_iter[i] = tos; + topo_next[i] = *tos; tos = i_iter; } } - if (i != r) - layer.push_or_transit_list.erase(layer_list_iter[i]); // ?? - else + if (i != r) { + i = topo_next[i]; + ++current[i]; + } else break; } } // while (1) - } // if - } // for + } // if color == white ... - if (bos != NULL) { + i_iter = vertices(g).first + i; + } // for all vertices in g + + if (bos != vertices(g).second) { i_iter = tos; do { i = *i_iter; ai = out_edges(i, g).first; a = *ai; - while (excess_flow[i] > 0) { - if (capacity[a] == 0 && residual_capacity[a] > 0) { - delta = std::min(excess_flow[i], residual_capacity[a]); - residual_capacity[a] -= delta; - residual_capacity[sister[a]] += delta; - excess_flow[i] -= delta; - excess_flow[target(a, g)] += delta; - } + while (excess_flow[i] > 0 && ai != out_edges(i, g).second) { + if (capacity[a] == 0 && residual_capacity[a] > 0) + push(a); ++ai; - } // while (excess_flow[i] > 0) - if (i == bos) + } + if (i == *bos) break; else - ++i_iter; + i = topo_next[i]; } while ( 1 ); } - + } // convert_preflow_to_flow() //======================================================================= @@ -368,9 +451,10 @@ namespace boost { // need to use random_access_property_map with these std::vector< FlowValue > excess_flow; - std::vector< typename Layer::list_iterator > layer_list_iter; + std::vector< list_node* > layer_list_ptr; std::vector< out_edge_iterator > current; std::vector< rank_size_type > rank; + std::vector color; // Edge Property Maps that must be interior to the graph SisterMap sister; @@ -382,7 +466,10 @@ namespace boost { rank_size_type lmin_push; // minimal layer with excess node std::queue Q; - enum { global_update_frequency = 1 }; + inline double global_update_frequency() { return 1.0; } + + list_node_next next; + list_node_prev prev; }; } // namespace detail @@ -396,33 +483,39 @@ namespace boost { EdgeCapacityMap cap, VertexIndexMap index_map, FlowValue& flow) { - detail::preflow_push + detail::push_relabel algo(g, cap, src, sink, index_map); algo.global_rank_update(); - num_relabels = 0; + typename graph_traits::vertices_size_type + num_relabels = 0, n = num_vertices(g); while (algo.lmax_push >= algo.lmin_push) { Layer& layer = algo.layers[algo.lmax_push]; - typename List::list_iterator i_iter = layer.push_list.begin(); - if (i_iter == layer.push_list.end()) + list_node* i_node = next(layer.push_list); + if (i_node == layer.push_list) --algo.lmax_push; else { - layer.push_list.pop_front(); - if (algo.push(*i_iter)) { // i must be relabelled - algo.relabel(*i_iter); + dlist_remove(next(layer.push_list), next, prev); + if (algo.discharge(*i_node)) { // returns true if done with out-edges + algo.relabel(*i_node); ++num_relabels; - if (layer.push_list.empty() && layer.transit_list.empty()) + + // Gap relabelling heuristic + if (dlist_empty(layer.push_list) && dlist_empty(layer.transit_list)) algo.gap(layer); - - if (n_r > global_update_frequency * n) { + + // Global relabelling heuristic + if (num_relabels > algo.global_update_frequency() * n) { algo.global_rank_update(); num_relabels = 0; } + } else { - layer.transit_list.erase(i_iter); + dlist_remove(i_node, next, prev); + delete i_node; } } @@ -495,7 +588,7 @@ namespace boost { // Vertex Properties: // - // layer_list_iter[v] points into push_list or trans_list + // layer_list_ptr[v] points into push_list or trans_list // current[v] (value_type == out_edge_iterator) // rank[v] also known as label // excess_flow[v]