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]