Port graph histograms and misc statistics to new filtering engine

Python binding will follow

src/graph/stats/Makefile.am

+## Process this file with automake to produce Makefile.in
+
+AM_CPPFLAGS =\
+    -I. -I.. \
+    -I../boost-workaround \
+    -DHAVE_CONFIG_H
+
+AM_CFLAGS=$(AM_CXXFLAGS)
+
+libgraph_tool_statsdir = $(pythondir)/graph_tool/stats
+
+libgraph_tool_stats_LTLIBRARIES = libgraph_tool_stats.la
+
+libgraph_tool_stats_la_includedir = $(pythondir)/graph_tool/include
+
+libgraph_tool_stats_la_SOURCES = \
+    graph_components.cc \
+    graph_histograms.cc \
+    graph_parallel.cc \
+    graph_stats_bind.cc
+
+libgraph_tool_stats_la_include_HEADERS = \
+    graph_components.hh \
+    graph_parallel.hh \
+    graph_histograms.hh
+
+libgraph_tool_stats_la_LIBADD = \
+    $(PYTHON_LDFLAGS) \
+    $(BOOST_LDFLAGS) \
+    $(OPENMP_LDFLAGS) \
+    -lboost_python \
+    -lboost_iostreams \
+    -lexpat
+
+# needed for typeinfo objects to work across DSO boundaries.
+# see http://gcc.gnu.org/faq.html#dso
+libgraph_tool_stats_la_LDFLAGS = \
+    -module \
+    -avoid-version \
+    -export-dynamic \
+    -no-undefined \
+    -Wl,-E

src/graph/stats/graph_components.cc

+// Copyright (C) 2008  Tiago de Paula Peixoto <tiago@forked.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_components.hh"
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void do_label_components(GraphInterface& gi, const string& prop)
+{
+    typedef property_map_types::apply<scalar_types,
+                                      GraphInterface::vertex_index_map_t>::type
+        vertex_props;
+    run_action<>()(gi, label_components(), vertex_props())
+        (vertex_prop(prop, gi));
+}
+
+void export_components()
+{
+    python::def("label_components", &do_label_components);
+};

src/graph/stats/graph_components.hh

+// Copyright (C) 2008  Tiago de Paula Peixoto <tiago@forked.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/>.
+
+#ifndef GRAPH_COMPONENTS_HH
+#define GRAPH_COMPONENTS_HH
+
+#include <boost/graph/connected_components.hpp>
+#include <boost/graph/strong_components.hpp>
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+// this will label the components of a graph to a given vertex property, from
+// [0, number of components - 1]. If the graph is directed the "strong
+// components" are used.
+struct label_components
+{
+    template <class Graph, class CompMap>
+    void operator()(const Graph* gp, CompMap comp_map) const
+    {
+        const Graph& g = *gp;
+        typedef typename graph_traits<Graph>::directed_category
+            directed_category;
+        get_components(g, comp_map,
+                       typename is_convertible<directed_category,
+                                               directed_tag>::type());
+    }
+
+    template <class Graph, class CompMap>
+    void get_components(const Graph& g, CompMap comp_map,
+                        boost::true_type is_directed) const
+    {
+        strong_components(g, comp_map);
+    }
+
+    template <class Graph, class CompMap>
+    void get_components(const Graph& g, CompMap comp_map,
+                        boost::false_type is_directed) const
+    {
+        connected_components(g, comp_map);
+    }
+};
+
+} // graph_tool namespace
+
+#endif // GRAPH_COMPONENTS_HH

src/graph/stats/graph_histograms.cc

@@ -13,3 +13,60 @@
 // 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_histograms.hh"
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+// this will return the vertex histogram of degrees or scalar properties
+python::object
+get_vertex_histogram(const GraphInterface& gi, GraphInterface::deg_t deg,
+                     const vector<long double>& bins)
+{
+    python::object hist;
+    python::object ret_bins;
+
+    run_action<>()(gi, get_histogram<VertexHistogramFiller>(hist, bins,
+                                                            ret_bins),
+                   all_selectors())(degree_selector(deg, gi));
+    return python::make_tuple(hist, ret_bins);
+}
+
+// this will return the vertex histogram of degrees or scalar properties
+python::object
+get_edge_histogram(GraphInterface& gi, const string& prop,
+                   const vector<long double>& bins)
+{
+    python::object hist;
+    python::object ret_bins;
+
+    bool directed = gi.GetDirected();
+    gi.SetDirected(false);
+
+    typedef property_map_types::apply<scalar_types,
+                                      GraphInterface::edge_index_map_t,
+                                      mpl::bool_<true> >::type
+        edge_props;
+
+    run_action<graph_tool::detail::always_directed>()
+        (gi, get_histogram<EdgeHistogramFiller>(hist, bins, ret_bins),
+         edge_props())(edge_prop(prop,gi));
+    gi.SetDirected(directed);
+
+    return python::make_tuple(hist, ret_bins);
+}
+
+using namespace boost::python;
+
+void export_histograms()
+{
+    def("get_vertex_histogram", &get_vertex_histogram);
+    def("get_edge_histogram", &get_edge_histogram);
+}

src/graph/stats/graph_histograms.hh

@@ -13,223 +13,168 @@
 // You should have received a copy of the GNU General Public License
 // along with this program. If not, see <http://www.gnu.org/licenses/>.
 
-// generalized functor to obtain histogram of different types of "degrees"
-struct get_vertex_histogram
+#ifndef GRAPH_HISTOGRAMS_HH
+#define GRAPH_HISTOGRAMS_HH
+
+#include <algorithm>
+#include <boost/numeric/conversion/bounds.hpp>
+#include <boost/numeric/conversion/cast.hpp>
+#include <boost/python/object.hpp>
+#include <boost/python/list.hpp>
+#include <boost/python/extract.hpp>
+#include "numpy_bind.hh"
+#include "histogram.hh"
+#include "shared_map.hh"
+
+namespace graph_tool
 {
-    template <class Graph, class DegreeSelector, class Hist>
-    void operator()(const Graph* gp, DegreeSelector deg, Hist& hist) const
-    {
-        const Graph& g = *gp;
-        SharedMap<Hist> s_hist(hist);
-        int i, N = num_vertices(g);
-        #pragma omp parallel for default(shared) private(i) \
-            firstprivate(s_hist) 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;
-            s_hist[deg(v,g)]++;
-        }
-        s_hist.Gather();
-    }
-};
-
+using namespace std;
+using namespace boost;
 
-// this will return the vertex histogram of degrees or scalar properties
-hist_t
-GraphInterface::GetVertexHistogram(GraphInterface::deg_t deg) const
+class VertexHistogramFiller
 {
-    hist_t hist;
-
-    try
+public:
+    template <class Graph, class DegreeSelector, class ValueType>
+    void UpdateRange(Graph& g,
+                     typename graph_traits<Graph>::vertex_descriptor v,
+                     DegreeSelector& deg, pair<ValueType,ValueType>& range)
     {
-        run_action<>()(*this, bind<void>(get_vertex_histogram(), _1, _2,
-                                         var(hist)),
-                       all_selectors())(degree_selector(deg, _properties));
+        ValueType val = deg(v, g);
+        range.first = min(range.first, val);
+        range.second = max(range.second, val);
     }
-    catch (dynamic_get_failure &e)
+
+    template <class Graph, class DegreeSelector, class Hist>
+    void operator()(Graph& g, typename graph_traits<Graph>::vertex_descriptor v,
+                    DegreeSelector& deg, Hist& hist)
     {
-        throw GraphException("error getting scalar property: " +
-                             string(e.what()));
+        typename Hist::point_t p;
+        p[0] = deg(v, g);
+        hist.PutValue(p);
     }
+};
 
-    return hist;
-}
-
-// generalized functor to obtain histogram of edge properties
-struct get_edge_histogram
+class EdgeHistogramFiller
 {
-    template <class Graph, class Prop, class Hist>
-    void operator()(const Graph* gp, Prop eprop, Hist& hist) const
+public:
+    template <class Graph, class EdgeProperty, class ValueType>
+    void UpdateRange(Graph& g,
+                     typename graph_traits<Graph>::vertex_descriptor v,
+                     EdgeProperty& eprop, pair<ValueType,ValueType>& range)
     {
-        const Graph& g = *gp;
-        SharedMap<Hist> s_hist(hist);
-
-        int i, N = num_vertices(g);
-        #pragma omp parallel for default(shared) private(i) \
-            firstprivate(s_hist)  schedule(dynamic)
-        for (i = 0; i < N; ++i)
+        typename graph_traits<Graph>::out_edge_iterator e, e_begin, e_end;
+        tie(e_begin,e_end) = out_edges(v, g);
+        for(e = e_begin; e != e_end; ++e)
         {
-            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
-            if (v == graph_traits<Graph>::null_vertex())
-                continue;
-
-            typename graph_traits<Graph>::out_edge_iterator e, e_begin, e_end;
-            tie(e_begin,e_end) = out_edges(v,g);
-            for(e = e_begin; e != e_end; ++e)
-                s_hist[eprop[*e]]++;
+            ValueType val = eprop[*e];
+            range.first = min(range.first, val);
+            range.second = max(range.second, val);
         }
-        s_hist.Gather();
+    }
 
-        typedef typename graph_traits<Graph>::directed_category
-            directed_category;
-        if(is_convertible<directed_category,undirected_tag>::value)
+    template <class Graph, class EdgeProperty, class Hist>
+    void operator()(Graph& g, typename graph_traits<Graph>::vertex_descriptor v,
+                    EdgeProperty& eprop, Hist& hist)
+    {
+        typename graph_traits<Graph>::out_edge_iterator e, e_begin, e_end;
+        tie(e_begin,e_end) = out_edges(v,g);
+        for(e = e_begin; e != e_end; ++e)
         {
-            for (typeof(s_hist.begin()) iter = s_hist.begin();
-                 iter != s_hist.end(); ++iter)
-                iter->second /= typename Hist::value_type::second_type(2);
+            typename Hist::point_t p;
+            p[0] = eprop[*e];
+            hist.PutValue(p);
         }
     }
 };
 
-
-// returns the histogram of a given edge property
-hist_t GraphInterface::GetEdgeHistogram(string property) const
+// generalized functor to obtain histogram of different types of "degrees"
+template <class HistogramFiller>
+struct get_histogram
 {
-    hist_t hist;
-    try
-    {
-        run_action<>()(*this, bind<void>(get_edge_histogram(), _1, _2,
-                                         var(hist)),
-                       edge_scalar_properties())
-            (prop(property, _edge_index, _properties));
-    }
-    catch (dynamic_get_failure& e)
+    get_histogram(python::object& hist, const vector<long double>& bins,
+                  python::object& ret_bins)
+        : _hist(hist), _bins(bins), _ret_bins(ret_bins) {}
+
+    template <class Graph, class DegreeSelector>
+    void operator()(Graph* gp, DegreeSelector deg) const
     {
-        throw GraphException("error getting scalar property: " +
-                             string(e.what()));
-    }
+        Graph& g = *gp;
 
-    return hist;
-}
+        typedef typename DegreeSelector::value_type value_type;
+        typedef Histogram<value_type, size_t, 1> hist_t;
 
-// this will label the components of a graph to a given vertex property, from
-// [0, number of components - 1]. If the graph is directed the "strong
-// components" are used.
-struct label_components
-{
-    template <class Graph, class CompMap>
-    void operator()(const Graph* gp, CompMap comp_map) const
-    {
-        const Graph& g = *gp;
-        typedef typename graph_traits<Graph>::directed_category
-            directed_category;
-        get_components(g, comp_map,
-                       typename is_convertible<directed_category,
-                                               directed_tag>::type());
-    }
+        HistogramFiller filler;
 
-    template <class Graph, class CompMap>
-    void get_components(const Graph& g, CompMap comp_map,
-                        boost::true_type is_directed) const
-    {
-        strong_components(g, comp_map);
-    }
+        vector<value_type> bins;
+        for (size_t i = 0; i < bins.size(); ++i)
+        {
+            // we'll attempt to recover from out of bounds conditions
+            try
+            {
+                bins[i] = numeric_cast<value_type,long double>(_bins[i]);
+            }
+            catch (boost::numeric::negative_overflow&)
+            {
+                bins[i] = boost::numeric::bounds<value_type>::lowest();
+            }
+            catch (boost::numeric::positive_overflow&)
+            {
+                bins[i] = boost::numeric::bounds<value_type>::highest();
+            }
+        }
 
-    template <class Graph, class CompMap>
-    void get_components(const Graph& g, CompMap comp_map,
-                        boost::false_type is_directed) const
-    {
-        connected_components(g, comp_map);
-    }
+        // sort the bins
+        sort(bins.begin(), bins.end());
 
-};
+        // clean bins of zero size
+        vector<value_type> temp_bin(1);
+        temp_bin[0] = bins[0];
+        for (size_t j = 1; j < bins.size(); ++j)
+        {
+            if (bins[j] > bins[j-1])
+                temp_bin.push_back(bins[j]);
+        }
+        bins = temp_bin;
 
-void GraphInterface::LabelComponents(string prop)
-{
-    try
-    {
-        run_action<>()(*this, label_components(), _1, vertex_scalar_selectors())
-            (degree_selector(prop, _properties));
-    }
-    catch (property_not_found)
-    {
-        typedef vector_property_map<int32_t, vertex_index_map_t> comp_map_t;
-        comp_map_t comp_map(_vertex_index);
-        _properties.property(prop, comp_map);
-        run_action<>()(*this, bind<void>(label_components(), _1, comp_map))();
-    }
-}
+        // find the data range
+        pair<value_type,value_type> range;
 
-// label parallel edges in the order they are found, starting from 0
-struct label_parallel_edges
-{
-    template <class Graph, class EdgeIndexMap, class ParallelMap>
-    void operator()(const Graph* gp, EdgeIndexMap edge_index,
-                    ParallelMap parallel) const
-    {
-        const Graph& g = *gp;
-        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
+        typename graph_traits<Graph>::vertex_iterator vi,vi_end;
+        range.first = boost::numeric::bounds<value_type>::highest();
+        range.second = boost::numeric::bounds<value_type>::lowest();
+        for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
+            filler.UpdateRange(g, *vi, deg, range);
+
+        boost::array<pair<value_type, value_type>, 1> data_range;
+        data_range[0] = range;
+        boost::array<vector<value_type>, 1> bin_list;
+        bin_list[0] = bins;
+
+        hist_t hist(bin_list, data_range);
+        SharedHistogram<hist_t> s_hist(hist);
 
         int i, N = num_vertices(g);
-        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+        #pragma omp parallel for default(shared) private(i) \
+            firstprivate(s_hist) 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;
-
-            tr1::unordered_set<edge_t,DescriptorHash<EdgeIndexMap> >
-                p_edges(0, DescriptorHash<EdgeIndexMap>(edge_index));
-
-            typename graph_traits<Graph>::out_edge_iterator e1, e2, e_end;
-            for (tie(e1, e_end) = out_edges(v, g); e1 != e_end; ++e1)
-            {
-                if (p_edges.find(*e1) != p_edges.end())
-                    continue;
-                size_t n = 0;
-                put(parallel, *e1, n);
-                for (tie(e2, e_end) = out_edges(v, g); e2 != e_end; ++e2)
-                    if (*e2 != *e1 && target(*e1, g) == target(*e2, g))
-                    {
-                        put(parallel, *e2, ++n);
-                        p_edges.insert(*e2);
-                    }
-            }
+            filler(g, *vi, deg, s_hist);
         }
-    }
-};
+        s_hist.Gather();
 
-void GraphInterface::LabelParallelEdges(string property)
-{
-    try
-    {
-        run_action<>()(*this, bind<void>(label_parallel_edges(), _1,
-                                         _edge_index, _2),
-                       edge_scalar_properties())
-            (prop(property, _edge_index, _properties));
-    }
-    catch (property_not_found)
-    {
-        typedef vector_property_map<int32_t,edge_index_map_t> parallel_map_t;
-        parallel_map_t parallel_map(_edge_index);
-        _properties.property(property, parallel_map);
-        run_action<>()(*this, bind<void>(label_parallel_edges(), _1,
-                                         _edge_index, parallel_map))();
+        bin_list = hist.GetBins();
+        python::object ret_bins = wrap_vector_owned(bin_list[0]);
+        _ret_bins = ret_bins;
+        _hist = wrap_multi_array_owned<size_t,1>(hist.GetArray());
     }
-}
-
-
-// this inserts the edge index map as a property
-void GraphInterface::InsertEdgeIndexProperty(string property)
-{
-    _properties.property(property, _edge_index);
-}
+    python::object& _hist;
+    const vector<long double>& _bins;
+    python::object& _ret_bins;
+};
 
-// this inserts the vertex index map as a property
-void GraphInterface::InsertVertexIndexProperty(string property)
-{
-    _properties.property(property, _vertex_index);
-}
+} // graph_tool namespace
 
+#endif // GRAPH_HISTOGRAMS_HH

src/graph/stats/graph_parallel.cc

+// Copyright (C) 2008  Tiago de Paula Peixoto <tiago@forked.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_parallel.hh"
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void do_label_parallel_edges(GraphInterface& gi, const string& property)
+{
+    run_action<>()(gi, bind<void>(label_parallel_edges(), _1,
+                                  gi.GetEdgeIndex(), _2),
+                   edge_scalar_properties())
+            (edge_prop(property, gi));
+}