Commit b7552219 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Implement correct behavior for min_st_cut()

This fixes ticket #167
parent 31373d8f
......@@ -27,51 +27,45 @@ void augment_graph(Graph& g, AugmentedMap augmented, CapacityMap capacity,
ReversedMap rmap, ResidualMap res,
bool detect_reversed = false)
{
typename graph_traits<Graph>::edge_iterator e, e_end;
for (tie(e,e_end) = edges(g); e != e_end; ++e)
augmented[*e] = 0;
for (auto e : edges_range(g))
augmented[e] = 0;
typename graph_traits<Graph>::vertex_iterator v, v_end;
vector<typename graph_traits<Graph>::edge_descriptor> e_list;
for (tie(v,v_end) = vertices(g); v != v_end; ++v)
for (auto v : vertices_range(g))
{
e_list.clear();
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e, e_end) = out_edges(*v, g); e != e_end; ++e)
for (auto e : out_edges_range(v, g))
{
if (detect_reversed && augmented[*e] == 0)
if (detect_reversed && augmented[e] == 0)
{
typename graph_traits<Graph>::out_edge_iterator e2, e2_end;
for (tie(e2, e2_end) = out_edges(target(*e, g), g);
e2 != e2_end; ++e2)
for (auto e2 : out_edges_range(target(e, g), g))
{
if (augmented[*e2] != 0)
if (augmented[e2] != 0)
continue;
if (target(*e2, g) == *v)
if (target(e2, g) == v)
{
augmented[*e] = 2;
augmented[*e2] = 2;
rmap[*e] = *e2;
rmap[*e2] = *e;
augmented[e] = 2;
augmented[e2] = 2;
rmap[e] = e2;
rmap[e2] = e;
break;
}
}
}
if (augmented[*e] == 0)
e_list.push_back(*e);
if (augmented[e] == 0)
e_list.push_back(e);
}
for (size_t i = 0; i < e_list.size(); ++i)
for (auto& e : e_list)
{
typename graph_traits<Graph>::edge_descriptor ae;
ae = add_edge(target(e_list[i],g), source(e_list[i],g), g).first;
auto ae = add_edge(target(e, g), source(e, g), g).first;
augmented[ae] = 1;
capacity[ae] = 0;
rmap[e_list[i]] = ae;
rmap[ae] = e_list[i];
rmap[e] = ae;
rmap[ae] = e;
res[ae] = 0;
}
}
......@@ -81,19 +75,37 @@ template <class Graph, class AugmentedMap>
void deaugment_graph(Graph& g, AugmentedMap augmented)
{
vector<typename graph_traits<Graph>::edge_descriptor> e_list;
typename graph_traits<Graph>::vertex_iterator v, v_end;
for (tie(v,v_end) = vertices(g); v != v_end; ++v)
for (auto v : vertices_range(g))
{
e_list.clear();
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e, e_end) = out_edges(*v, g); e != e_end; ++e)
for (auto e : out_edges_range(v, g))
{
if (augmented[*e] == 1)
e_list.push_back(*e);
if (augmented[e] == 1)
e_list.push_back(e);
}
for (size_t i = 0; i < e_list.size(); ++i)
remove_edge(e_list[i], g);
for (auto& e : e_list)
remove_edge(e, g);
}
}
template <class Graph, class CapacityMap, class ResidualMap,
class AugmentedMap>
void residual_graph(Graph& g, CapacityMap capacity, ResidualMap res,
AugmentedMap augmented)
{
vector<typename graph_traits<Graph>::edge_descriptor> e_list;
for (auto e : edges_range(g))
{
if (capacity[e] - res[e] > 0)
e_list.push_back(e);
}
for (auto& e : e_list)
{
auto ne = add_edge(target(e, g), source(e, g), g);
augmented[ne.first] = true;
}
}
......
......@@ -26,6 +26,8 @@ void kolmogorov_max_flow(GraphInterface& gi, size_t src, size_t sink,
boost::any capacity, boost::any res);
bool max_cardinality_matching(GraphInterface& gi, boost::any match);
double min_cut(GraphInterface& gi, boost::any weight, boost::any part_map);
void get_residual_graph(GraphInterface& gi, boost::any capacity, boost::any res,
boost::any oaugment);
#include <boost/python.hpp>
using namespace boost::python;
......@@ -37,4 +39,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_flow)
def("kolmogorov_max_flow", &kolmogorov_max_flow);
def("max_cardinality_matching", &max_cardinality_matching);
def("min_cut", &min_cut);
def("residual_graph", &get_residual_graph);
}
......@@ -19,6 +19,7 @@
#include "graph.hh"
#include "graph_properties.hh"
#include "graph_augment.hh"
#include <boost/graph/stoer_wagner_min_cut.hpp>
using namespace std;
......@@ -60,3 +61,27 @@ double min_cut(GraphInterface& gi, boost::any weight, boost::any part_map)
weight_maps(), writable_vertex_scalar_properties())(weight, part_map);
return mc;
}
struct do_get_residual_graph
{
template <class Graph, class CapacityMap, class ResidualMap,
class AugmentedMap>
void operator()(Graph& g, CapacityMap capacity, ResidualMap res,
AugmentedMap augmented) const
{
residual_graph(g, capacity, res, augmented);
}
};
void get_residual_graph(GraphInterface& gi, boost::any capacity,
boost::any res, boost::any oaugment)
{
typedef property_map_type::apply<uint8_t,
GraphInterface::edge_index_map_t>::type
emap_t;
emap_t augment = boost::any_cast<emap_t>(oaugment);
run_action<>()
(gi, std::bind(do_get_residual_graph(), placeholders::_1,
placeholders::_2, placeholders::_3, augment),
edge_scalar_properties(), edge_scalar_properties())(capacity, res);
}
......@@ -347,7 +347,7 @@ def boykov_kolmogorov_max_flow(g, source, target, capacity, residual=None):
_prop("e", g, residual))
return residual
def min_st_cut(g, source, residual):
def min_st_cut(g, source, capacity, residual):
r"""
Get the minimum source-target cut, given the residual capacity of the edges.
......@@ -357,13 +357,13 @@ def min_st_cut(g, source, residual):
Graph to be used.
source : Vertex
The source vertex.
capacity : :class:`~graph_tool.PropertyMap`
Edge property map with the edge capacities.
residual : :class:`~graph_tool.PropertyMap`
Edge property map where the residual capacity is stored.
Edge property map with the residual capacities (capacity - flow).
Returns
-------
min_cut : float
The value of the minimum cut.
partition : :class:`~graph_tool.PropertyMap`
Boolean-valued vertex property map with the cut partition. Vertices with
value `True` belong to the source side of the cut.
......@@ -372,38 +372,42 @@ def min_st_cut(g, source, residual):
-----
The source-side of the cut set is obtained by following all vertices which
are reachable from the source via edges with nonzero residual capacity.
are reachable from the source in the residual graph (i.e. via edges
with nonzero residual capacity, and reversed edges with nonzero flow).
This algorithm runs in :math:`O(V+E)` time.
Examples
--------
>>> g = gt.load_graph("flow-example.xml.gz")
>>> cap = g.edge_properties["cap"]
>>> src, tgt = g.vertex(0), g.vertex(1)
>>> g = gt.load_graph("mincut-st-example.xml.gz")
>>> cap = g.edge_properties["weight"]
>>> src, tgt = g.vertex(0), g.vertex(7)
>>> res = gt.boykov_kolmogorov_max_flow(g, src, tgt, cap)
>>> mc, part = gt.min_st_cut(g, src, res)
>>> part = gt.min_st_cut(g, src, cap, res)
>>> mc = sum([cap[e] - res[e] for e in g.edges() if part[e.source()] != part[e.target()]])
>>> print(mc)
14.331937627198545
3
>>> pos = g.vertex_properties["pos"]
>>> res.a = cap.a - res.a # the actual flow
>>> gt.graph_draw(g, pos=pos, edge_pen_width=gt.prop_to_size(res, mi=0, ma=3, power=1),
... vertex_fill_color=part, output="example-min-st-cut.pdf")
>>> gt.graph_draw(g, pos=pos, edge_pen_width=gt.prop_to_size(cap, mi=3, ma=10, power=1),
... edge_text=res, vertex_fill_color=part, vertex_text=g.vertex_index,
... vertex_font_size=18, edge_font_size=18, output="example-min-st-cut.pdf")
<...>
.. testcode::
:hide:
gt.graph_draw(g, pos=pos, edge_pen_width=gt.prop_to_size(res, mi=0, ma=3, power=1), vertex_fill_color=part,
output="example-min-st-cut.png")
gt.graph_draw(g, pos=pos, edge_pen_width=gt.prop_to_size(cap, mi=3, ma=10, power=1),
edge_text=res, vertex_fill_color=part, vertex_text=g.vertex_index,
vertex_font_size=18, edge_font_size=18, output="example-min-st-cut.png")
.. figure:: example-min-st-cut.*
:align: center
Edge flows obtained with the Boykov-Kolmogorov algorithm. The source and
target are on the lower left and upper right corners, respectively. The
edge flows are represented by the edge width. Vertices of the same color
are on the same side of a minimum cut.
target are labeled ``0`` and ``7``, respectively. The edge capacities are
represented by the edge width, and the maximum flow by the edge labels.
Vertices of the same color are on the same side the minimum cut.
References
----------
......@@ -411,13 +415,22 @@ def min_st_cut(g, source, residual):
"""
if not g.is_directed():
raise ValueError("The graph provided must be directed!")
augment = g.new_edge_property("bool")
libgraph_tool_flow.residual_graph(g._Graph__graph,
_prop("e", g, capacity),
_prop("e", g, residual),
_prop("e", g, augment))
em = g.new_edge_property("bool")
em.a = residual.a[:len(em.a)] > 0
em.a = (residual.a[:len(em.a)] > 0) + augment.a
u = GraphView(g, efilt=em)
part = label_out_component(u, source)
g.own_property(part)
mc = sum(residual[e] for e in source.out_edges())
return mc, part
part = g.own_property(part)
# cleanup augmented edges
u = GraphView(g, efilt=augment)
u.clear_edges()
return part
def min_cut(g, weight):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment