Implement pseudo_diameter()

src/graph/topology/Makefile.am

@@ -18,6 +18,7 @@ libgraph_tool_topology_la_SOURCES = \
     graph_all_distances.cc \
     graph_components.cc \
     graph_distance.cc \
+    graph_diameter.cc \
     graph_dominator_tree.cc \
     graph_isomorphism.cc \
     graph_minimum_spanning_tree.cc \

src/graph/topology/graph_diameter.cc

+// Copyright (C) 2006-2011 Tiago de Paula Peixoto <tiago@skewed.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_properties.hh"
+#include "graph_selectors.hh"
+
+#include <boost/graph/breadth_first_search.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+template <class DistMap>
+class bfs_diam_visitor:
+    public boost::bfs_visitor<null_visitor>
+{
+public:
+    bfs_diam_visitor(DistMap dist_map, size_t& v)
+        : _dist_map(dist_map), _v(v), _dist(0), _max_dist(0),
+          _min_k(numeric_limits<size_t>::max()) {}
+
+    template <class Graph>
+    void tree_edge(typename graph_traits<Graph>::edge_descriptor e,
+                   Graph& g)
+    {
+        typename graph_traits<Graph>::vertex_descriptor v = target(e,g);
+        size_t dist = _dist_map[source(e,g)] + 1;
+        if ((dist > _max_dist) ||
+            (dist == _max_dist && total_degreeS()(v, g) <= _min_k))
+        {
+            _max_dist = dist;
+            _min_k = total_degreeS()(v, g);
+            _v = v;
+        }
+        _dist_map[v] = dist;
+    }
+
+private:
+    DistMap _dist_map;
+    size_t& _v;
+    size_t _dist;
+    size_t _max_dist;
+    size_t _min_k;
+};
+
+template <class DistMap>
+class djk_diam_visitor:
+    public boost::dijkstra_visitor<null_visitor>
+{
+public:
+    djk_diam_visitor(DistMap dist_map, size_t& v)
+        : _dist_map(dist_map), _v(v), _max_dist(0),
+          _min_k(numeric_limits<size_t>::max()) {}
+
+    template <class Graph>
+    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor v,
+                        Graph& g)
+    {
+        if ((_dist_map[v] > _max_dist) ||
+            (_dist_map[v] == _max_dist && total_degreeS()(v, g) <= _min_k))
+        {
+            _max_dist = _dist_map[v];
+            _min_k = total_degreeS()(v, g);
+            _v = v;
+        }
+    }
+
+private:
+    DistMap _dist_map;
+    size_t& _v;
+    typename property_traits<DistMap>::value_type _max_dist;
+    size_t _min_k;
+
+
+};
+
+
+struct do_bfs_search
+{
+    template <class Graph, class VertexIndexMap>
+    void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
+                    size_t& target, long double& max_dist) const
+    {
+        typedef unchecked_vector_property_map<size_t, VertexIndexMap> dist_map_t;
+        dist_map_t dist_map(vertex_index, num_vertices(g));
+        int i, N = num_vertices(g);
+        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+        for (i = 0; i < N; ++i)
+        {
+            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
+            if (v == graph_traits<Graph>::null_vertex())
+                continue;
+            dist_map[v] = numeric_limits<size_t>::max();
+        }
+        dist_map[vertex(source,g)] = 0;
+
+        unchecked_vector_property_map<boost::default_color_type, VertexIndexMap>
+            color_map(vertex_index, num_vertices(g));
+        target = source;
+        breadth_first_search(g, vertex(source, g),
+                             visitor(bfs_diam_visitor<dist_map_t>
+                                     (dist_map, ref(target))).
+                             vertex_index_map(vertex_index).
+                             color_map(color_map));
+        max_dist = dist_map[vertex(target, g)];
+    }
+};
+
+struct do_djk_search
+{
+    template <class Graph, class VertexIndexMap, class WeightMap>
+    void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
+                    WeightMap weight, size_t& target, long double& max_dist) const
+    {
+        typedef unchecked_vector_property_map<typename property_traits<WeightMap>::value_type,
+                                              VertexIndexMap> dist_map_t;
+        dist_map_t dist_map(vertex_index, num_vertices(g));
+        target = source;
+        dijkstra_shortest_paths(g, vertex(source, g),
+                                weight_map(weight).
+                                distance_map(dist_map).
+                                vertex_index_map(vertex_index).
+                                visitor(djk_diam_visitor<dist_map_t>
+                                        (dist_map, ref(target))));
+        max_dist = dist_map[vertex(target, g)];
+    }
+};
+
+python::object get_diam(GraphInterface& gi, size_t source, boost::any weight)
+{
+    size_t target;
+    long double max_dist;
+    if (weight.empty())
+    {
+        run_action<>()
+            (gi, bind<void>(do_bfs_search(), _1, source, gi.GetVertexIndex(),
+                            ref(target), ref(max_dist)))();
+    }
+    else
+    {
+        run_action<>()
+            (gi, bind<void>(do_djk_search(), _1, source, gi.GetVertexIndex(),
+                            _2, ref(target), ref(max_dist)),
+             edge_scalar_properties())(weight);
+
+    }
+    return python::make_tuple(target, max_dist);
+}
+
+void export_diam()
+{
+    python::def("get_diam", &get_diam);
+};

src/graph/topology/graph_topology.cc

@@ -42,6 +42,7 @@ void export_components();
 void export_similarity();
 void export_dists();
 void export_all_dists();
+void export_diam();
 
 BOOST_PYTHON_MODULE(libgraph_tool_topology)
 {
@@ -57,4 +58,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_topology)
     export_similarity();
     export_dists();
     export_all_dists();
+    export_diam();
 }

src/graph_tool/topology/__init__.py

@@ -30,6 +30,7 @@ Summary
 
    shortest_distance
    shortest_path
+   pseudo_diameter
    similarity
    isomorphism
    subgraph_isomorphism
@@ -58,7 +59,7 @@ __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "min_spanning_tree", "dominator_tree", "topological_sort",
            "transitive_closure", "label_components", "label_largest_component",
            "label_biconnected_components", "shortest_distance",
-           "shortest_path", "is_planar", "similarity"]
+           "shortest_path", "pseudo_diameter", "is_planar", "similarity"]
 
 
 def similarity(g1, g2, label1=None, label2=None, norm=True):
@@ -844,8 +845,7 @@ def shortest_path(g, source, target, weights=None, pred_map=None):
     target : :class:`~graph_tool.Vertex`
         Target vertex of the search.
     weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
-        The edge weights. If provided, the minimum spanning tree will minimize
-        the edge weights.
+        The edge weights.
     pred_map :  :class:`~graph_tool.PropertyMap` (optional, default: None)
         Vertex property map with the predecessors in the search tree. If this is
         provided, the shortest paths are not computed, and are obtained directly
@@ -926,6 +926,76 @@ def shortest_path(g, source, target, weights=None, pred_map=None):
     return vlist, elist
 
 
+def pseudo_diameter(g, source=None, weights=None):
+    """
+    Compute the pseudo-diameter of the graph.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    source : :class:`~graph_tool.Vertex` (optional, default: `None`)
+        Source vertex of the search. If not supplied, the first vertex
+        in the graph will be chosen.
+    weights : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
+        The edge weights.
+
+    Returns
+    -------
+    pseudo_diameter : int
+        The pseudo-diameter of the graph.
+    end_points : pair of :class:`~graph_tool.Vertex`
+        The two vertices which correspond to the pseudo-diameter found.
+
+    Notes
+    -----
+
+    The pseudo-diameter is an approximate graph diameter. It is obtained by
+    starting from a vertex `source`, and finds a vertex `target` that is
+    farthest away from `source`. This process is repeated by treating
+    `target` as the new starting vertex, and ends when the graph distance no
+    longer increases. A vertex from the last level set that has the smallest
+    degree is chosen as the final starting vertex u, and a traversal is done
+    to see if the graph distance can be increased. This graph distance is
+    taken to be the pseudo-diameter.
+
+    The paths are computed with a breadth-first search (BFS) or Dijkstra's
+    algorithm [dijkstra]_, if weights are given.
+
+    The algorithm runs in :math:`O(V + E)` time, or :math:`O(V \log V)` if
+    weights are given.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, poisson
+    >>> seed(42)
+    >>> g = gt.random_graph(300, lambda: (poisson(3), poisson(3)))
+    >>> dist, ends = gt.pseudo_diameter(g)
+    >>> print dist
+    >>> print end
+
+    References
+    ----------
+    .. [pseudo-diameter] http://en.wikipedia.org/wiki/Distance_%28graph_theory%29
+    """
+
+    if source is None:
+        source = g.vertex(0)
+    dist, target = 0, source
+    while True:
+        new_source = target
+        new_target, new_dist = libgraph_tool_topology.get_diam(g._Graph__graph,
+                                                               int(new_source),
+                                                               _prop("e", g, weights))
+        if new_dist > dist:
+            target = new_target
+            source = new_source
+            dist = new_dist
+        else:
+            break
+    return dist, (g.vertex(source), g.vertex(target))
+
+
 def is_planar(g, embedding=False, kuratowski=False):
     """
     Test if the graph is planar.