Added (undocumented) functions for loop-erased random walk and random generation of spanning trees

[SVN r63334]

include/boost/graph/loop_erased_random_walk.hpp

@@ -0,0 +1,110 @@
+// Copyright 2010 The Trustees of Indiana University.
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+//  Authors: Jeremiah Willcock
+//           Andrew Lumsdaine
+
+#ifndef BOOST_GRAPH_LOOP_ERASED_RANDOM_WALK_HPP
+#define BOOST_GRAPH_LOOP_ERASED_RANDOM_WALK_HPP
+
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/graph/random.hpp>
+#include <boost/next_prior.hpp>
+#include <vector>
+#include <cassert>
+
+namespace boost {
+
+  // Do a loop-erased random walk from vertex s to any vertex colored black (or
+  // actually any color other than white or gray) in the color map.  The color
+  // white is for vertices that are not part of the path, while gray is for
+  // those that are on the path (for cycle detection).  The vector path is used
+  // for temporary storage and as the result of the algorithm; while all
+  // elements of the path except the last have their colors set to gray upon
+  // return.  Vertex s must start off colored white.
+  //
+  // Useful references:
+  // http://everything2.com/title/loop-erased+random+walk
+  // Wikipedia page on "Loop-Erased Random Walk"
+
+  template <typename Graph, typename ColorMap, typename NextEdge>
+  void loop_erased_random_walk(
+         const Graph& g,
+         typename boost::graph_traits<Graph>::vertex_descriptor s,
+         NextEdge next_edge,
+         ColorMap color,
+         std::vector<typename boost::graph_traits<Graph>::vertex_descriptor>& path
+  ) {
+    typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+    typedef typename boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
+    typedef typename boost::property_traits<ColorMap>::value_type color_t;
+    typedef boost::color_traits<color_t> color_gen;
+    
+    assert (get(color, s) == color_gen::white());
+    path.clear();
+    path.push_back(s);
+    put(color, s, color_gen::gray());
+    while (true) {
+      edge_descriptor e = next_edge(s, g);
+      vertex_descriptor t = target(e, g);
+      color_t t_color = get(color, t);
+      if (t_color == color_gen::white()) {
+        path.push_back(t);
+        put(color, t, color_gen::gray());
+        s = t;
+      } else if (t_color == color_gen::gray()) {
+        // Found a loop; delete from path from the first occurrence of t to the
+        // end, coloring vertices white.
+        typename std::vector<vertex_descriptor>::iterator it = std::find(path.begin(), path.end(), t);
+        assert (it != path.end());
+        ++it;
+        for (typename std::vector<vertex_descriptor>::iterator j = it; j != path.end(); ++j) {
+          put(color, *j, color_gen::white());
+        }
+        path.erase(it, path.end());
+        s = t;
+      } else {
+        // Done
+        path.push_back(t);
+        break;
+      }
+    }
+  }
+
+  template <typename Graph, typename Gen>
+  class unweighted_random_out_edge_gen {
+    Gen& gen;
+
+    typedef boost::graph_traits<Graph> gt;
+
+    public:
+    unweighted_random_out_edge_gen(Gen& gen): gen(gen) {}
+
+    typename gt::edge_descriptor
+    operator()(typename gt::vertex_descriptor src, const Graph& g) const {
+      return boost::random_out_edge(g, src, gen);
+    }
+  };
+
+  template <typename Graph, typename WeightMap, typename Gen>
+  class weighted_random_out_edge_gen {
+    WeightMap weight;
+    Gen& gen;
+
+    typedef boost::graph_traits<Graph> gt;
+
+    public:
+    weighted_random_out_edge_gen(const WeightMap& weight, Gen& gen): weight(weight), gen(gen) {}
+
+    typename gt::edge_descriptor
+    operator()(typename gt::vertex_descriptor src, const Graph& g) const {
+      return boost::weighted_random_out_edge(g, src, weight, gen);
+    }
+  };
+}
+
+#endif // BOOST_GRAPH_LOOP_ERASED_RANDOM_WALK_HPP

include/boost/graph/random.hpp

@@ -12,10 +12,12 @@
 
 #include <boost/graph/graph_traits.hpp>
 #include <boost/random/uniform_int.hpp>
+#include <boost/random/uniform_real.hpp>
 #include <boost/random/variate_generator.hpp>
 
 #include <boost/pending/property.hpp>
 #include <boost/graph/properties.hpp>
+#include <boost/graph/iteration_macros.hpp>
 #include <boost/next_prior.hpp>
 
 #include <boost/graph/adjacency_list.hpp>
@@ -24,6 +26,7 @@
 #include <boost/type_traits/is_convertible.hpp>
 
 #include <iostream>
+#include <cassert>
 
 namespace boost {
 
@@ -67,6 +70,43 @@ namespace boost {
       return *edges(g).first;
   }
 
+  template <typename Graph, typename RandomNumGen>
+  typename graph_traits<Graph>::edge_descriptor
+  random_out_edge(Graph& g, typename graph_traits<Graph>::vertex_descriptor src, RandomNumGen& gen) {
+    typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
+    typedef boost::uniform_int<degree_size_type> ui_type;
+    ui_type ui(0, out_degree(src, g) - 1);
+    boost::variate_generator<RandomNumGen&, ui_type>
+      variate(gen, ui);
+    typename graph_traits<Graph>::out_edge_iterator it = out_edges(src, g).first;
+    std::advance(it, variate());
+    return *it;
+  }
+
+  template <typename Graph, typename WeightMap, typename RandomNumGen>
+  typename graph_traits<Graph>::edge_descriptor
+  weighted_random_out_edge(Graph& g, typename graph_traits<Graph>::vertex_descriptor src, WeightMap weight, RandomNumGen& gen) {
+    typedef graph_traits<Graph> gt;
+    typedef typename gt::vertex_descriptor vertex_descriptor;
+    typedef typename property_traits<WeightMap>::value_type weight_type;
+    weight_type weight_sum(0);
+    BGL_FORALL_OUTEDGES_T(src, e, g, Graph) {weight_sum += get(weight, e);}
+    typedef boost::uniform_real<> ur_type;
+    ur_type ur(0, weight_sum);
+    boost::variate_generator<RandomNumGen&, ur_type>
+      variate(gen, ur);
+    weight_type chosen_weight = variate();
+    BGL_FORALL_OUTEDGES_T(src, e, g, Graph) {
+      weight_type w = get(weight, e);
+      if (chosen_weight < w) {
+        return e;
+      } else {
+        chosen_weight -= w;
+      }
+    }
+    assert (false); // Should not get here
+  }
+
   namespace detail {
     class dummy_property_copier {
     public:
@@ -200,5 +240,6 @@ namespace boost {
   
 }
 
+#include <boost/graph/iteration_macros_undef.hpp>
 
 #endif

include/boost/graph/random_spanning_tree.hpp

@@ -0,0 +1,103 @@
+// Copyright 2010 The Trustees of Indiana University.
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+//  Authors: Jeremiah Willcock
+//           Andrew Lumsdaine
+
+#ifndef BOOST_GRAPH_RANDOM_SPANNING_TREE_HPP
+#define BOOST_GRAPH_RANDOM_SPANNING_TREE_HPP
+
+#include <vector>
+#include <boost/graph/loop_erased_random_walk.hpp>
+#include <boost/graph/random.hpp>
+#include <boost/graph/iteration_macros.hpp>
+#include <boost/property_map/property_map.hpp>
+
+namespace boost {
+
+  namespace detail {
+    // Use Wilson's algorithm (based on loop-free random walks) to generate a
+    // random spanning tree.  The distribution of edges used is controlled by
+    // the next_edge() function, so this version allows either weighted or
+    // unweighted selection of trees.
+    // Algorithm is from http://en.wikipedia.org/wiki/Uniform_spanning_tree
+    template <typename Graph, typename PredMap, typename ColorMap, typename NextEdge>
+    void random_spanning_tree_internal(const Graph& g, PredMap pred, ColorMap color, NextEdge next_edge) {
+      typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+      typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
+
+      assert (num_vertices(g) >= 2); // g must also be undirected (or symmetric) and connected
+
+      typedef color_traits<typename property_traits<ColorMap>::value_type> color_gen;
+      BGL_FORALL_VERTICES_T(v, g, Graph) put(color, v, color_gen::white());
+
+      std::vector<vertex_descriptor> path;
+
+      vertex_descriptor s = *vertices(g).first;
+      put(color, s, color_gen::black());
+      put(pred, s, graph_traits<Graph>::null_vertex());
+
+      BGL_FORALL_VERTICES_T(v, g, Graph) {
+        if (get(color, v) != color_gen::white()) continue;
+        loop_erased_random_walk(g, v, next_edge, color, path);
+        for (typename std::vector<vertex_descriptor>::const_reverse_iterator i = path.rbegin();
+             boost::next(i) != path.rend(); ++i) {
+          typename std::vector<vertex_descriptor>::const_reverse_iterator j = i;
+          ++j;
+          assert (get(color, *j) == color_gen::gray());
+          put(color, *j, color_gen::black());
+          put(pred, *j, *i);
+        }
+      }
+    }
+  }
+
+  // Compute a uniformly-distributed spanning tree on a graph.
+  template <typename Graph, typename PredMap, typename ColorMap, typename Gen>
+  void random_spanning_tree(const Graph& g, PredMap pred, Gen& gen, ColorMap color) {
+    unweighted_random_out_edge_gen<Graph, Gen> random_oe(gen);
+    detail::random_spanning_tree_internal(g, pred, color, random_oe);
+  }
+
+  // Compute a uniformly-distributed spanning tree on a graph.
+  template <typename Graph, typename PredMap, typename Gen>
+  void random_spanning_tree(const Graph& g, PredMap pred, Gen& gen) {
+    std::vector<default_color_type> color_data(num_vertices(g));
+    random_spanning_tree(g, pred, gen, make_iterator_property_map(color_data.begin(), get(vertex_index, g)));
+  }
+
+  // Compute a weight-distributed spanning tree on a graph.
+  // Weighting works according to:
+  // @article{Mosbah1999263,
+  //   title = "Non-uniform random spanning trees on weighted graphs",
+  //   journal = "Theoretical Computer Science",
+  //   volume = "218",
+  //   number = "2",
+  //   pages = "263--271",
+  //   year = "1999",
+  //   note = "",
+  //   issn = "0304-3975",
+  //   doi = "DOI: 10.1016/S0304-3975(98)00325-9",
+  //   url = "http://www.sciencedirect.com/science/article/B6V1G-3WSV1D9-P/2/06bea092e23163c4884844cde4a5e92c",
+  //   author = "M. Mosbah and N. Saheb"
+  // }
+  template <typename Graph, typename PredMap, typename WeightMap, typename ColorMap, typename Gen>
+  void weighted_random_spanning_tree(const Graph& g, PredMap pred, WeightMap weight, Gen& gen, ColorMap color) {
+    weighted_random_out_edge_gen<Graph, WeightMap, Gen> random_oe(weight, gen);
+    detail::random_spanning_tree_internal(g, pred, color, random_oe);
+  }
+
+  // Compute a weight-distributed spanning tree on a graph.
+  template <typename Graph, typename PredMap, typename WeightMap, typename Gen>
+  void weighted_random_spanning_tree(const Graph& g, PredMap pred, WeightMap weight, Gen& gen) {
+    std::vector<default_color_type> color_data(num_vertices(g));
+    weighted_random_spanning_tree(g, pred, weight, gen, make_iterator_property_map(color_data.begin(), get(vertex_index, g)));
+  }
+}
+
+#include <boost/graph/iteration_macros_undef.hpp>
+
+#endif // BOOST_GRAPH_RANDOM_SPANNING_TREE_HPP

test/Jamfile.v2

@@ -116,6 +116,7 @@ test-suite graph_test :
     [ compile grid_graph_cc.cpp ]
     [ run grid_graph_test.cpp ]
     [ run incremental_components_test.cpp ]
+    [ run random_spanning_tree_test.cpp ../build//boost_graph ]
     [ run graphml_test.cpp ../build//boost_graph : : "graphml_test.xml" ]
     ;
 

test/random_spanning_tree_test.cpp

@@ -0,0 +1,69 @@
+// Copyright 2010 The Trustees of Indiana University.
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+//  Authors: Jeremiah Willcock
+//           Andrew Lumsdaine
+
+#include <boost/graph/random_spanning_tree.hpp>
+#include <boost/graph/grid_graph.hpp>
+#include <boost/array.hpp>
+#include <boost/random.hpp>
+#include <boost/property_map/shared_array_property_map.hpp>
+#include <boost/property_map/dynamic_property_map.hpp>
+#include <boost/graph/graphviz.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/graph/property_maps/constant_property_map.hpp>
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <boost/graph/iteration_macros.hpp>
+
+using namespace boost;
+using namespace std;
+
+typedef grid_graph<2> graph_type;
+typedef graph_traits<graph_type> gt;
+
+template <typename Graph, typename PredMap, typename WeightMap>
+void write_spanning_tree(const Graph& g, PredMap pred, WeightMap weight, string filename) {
+  shared_array_property_map<string, typename property_map<Graph, edge_index_t>::const_type> edge_style(num_edges(g), get(edge_index, g));
+  shared_array_property_map<string, typename property_map<Graph, vertex_index_t>::const_type> vertex_pos(num_vertices(g), get(vertex_index, g));
+  BGL_FORALL_EDGES_T(e, g, Graph) {
+    put(edge_style, e, (get(pred, target(e, g)) == source(e, g) || get(pred, source(e, g)) == target(e, g)) ? "bold" : "dotted");
+  }
+  BGL_FORALL_VERTICES_T(v, g, Graph) {
+    put(vertex_pos, v, lexical_cast<string>(v[0]) + "," + lexical_cast<string>(v[1]));
+  }
+  dynamic_properties dp;
+  dp.property("style", edge_style);
+  dp.property("pos", vertex_pos);
+  dp.property("len", weight);
+  dp.property("node_id", get(vertex_index, g));
+  ofstream out(filename.c_str());
+  write_graphviz_dp(out, g, dp);
+}
+
+int main(int, char**) {
+
+  array<size_t, 2> sizes = {{ 5, 5 }};
+  graph_type g(sizes);
+
+  shared_array_property_map<gt::vertex_descriptor, property_map<graph_type, vertex_index_t>::const_type> pred(num_vertices(g), get(vertex_index, g));
+  shared_array_property_map<double, property_map<graph_type, edge_index_t>::const_type> weight(num_edges(g), get(edge_index, g));
+
+  BGL_FORALL_EDGES(e, g, graph_type) {put(weight, e, (1. + get(edge_index, g, e)) / num_edges(g));}
+
+  mt19937 gen;
+  random_spanning_tree(g, pred, gen);
+  // write_spanning_tree(g, pred, constant_property_map<gt::edge_descriptor, double>(1.), "unweight_random_st.dot");
+  random_spanning_tree(g, pred, gen);
+  // write_spanning_tree(g, pred, constant_property_map<gt::edge_descriptor, double>(1.), "unweight_random_st2.dot");
+
+  weighted_random_spanning_tree(g, pred, weight, gen);
+  // write_spanning_tree(g, pred, weight, "weight_random_st.dot");
+
+  return 0;
+}