Enable vector types and python::object in vertex/edge_average()

src/graph/stats/graph_average.cc

@@ -23,6 +23,8 @@
 
 #include <boost/python.hpp>
 
+#include <boost/mpl/joint_view.hpp>
+
 using namespace std;
 using namespace boost;
 using namespace graph_tool;
@@ -31,25 +33,53 @@ using namespace graph_tool;
 python::object
 get_vertex_average(GraphInterface& gi, GraphInterface::deg_t deg)
 {
-    long double a, dev;
-    run_action<>()(gi, get_average<VertexAverageTraverse>(a,dev),
-                   scalar_selectors())(degree_selector(deg));
-    return python::make_tuple(a,dev);
+    python::object a, dev;
+    size_t count;
+
+    typedef mpl::joint_view<vertex_scalar_properties,
+                            vertex_scalar_vector_properties>
+        vertex_numeric_properties_c;
+
+    typedef property_map_types::apply<mpl::vector<python::object>,
+                                      GraphInterface::vertex_index_map_t>::type
+        python_properties;
+
+    typedef mpl::joint_view<vertex_numeric_properties_c,
+                            python_properties>
+        vertex_numeric_properties;
+
+    typedef mpl::transform<vertex_numeric_properties,
+                           scalar_selector_type,
+                           mpl::back_inserter<degree_selectors> >::type
+        numeric_selectors;
+
+    run_action<>()(gi, get_average<VertexAverageTraverse>(a,dev,count),
+                   numeric_selectors())(degree_selector(deg));
+    return python::make_tuple(a, dev, count);
 }
 
 // this will return the edge average of scalar properties
 python::object
 get_edge_average(GraphInterface& gi, boost::any prop)
 {
-    long double a, dev;
-    bool directed = gi.GetDirected();
-    gi.SetDirected(true);
-    run_action<graph_tool::detail::always_directed>()
-        (gi, get_average<EdgeAverageTraverse>(a, dev),
-         edge_scalar_properties())(prop);
-    gi.SetDirected(directed);
+    typedef mpl::joint_view<edge_scalar_properties,
+                            edge_scalar_vector_properties>
+        edge_numeric_properties_c;
+
+    typedef property_map_types::apply<mpl::vector<python::object>,
+                                      GraphInterface::edge_index_map_t>::type
+        python_properties;
 
-    return python::make_tuple(a, dev);
+    typedef mpl::joint_view<edge_numeric_properties_c,
+                            python_properties>
+        edge_numeric_properties;
+
+    python::object a, dev;
+    size_t count;
+    run_action<graph_tool::detail::always_directed>()
+        (gi, get_average<EdgeAverageTraverse>(a, dev, count),
+         edge_numeric_properties())(prop);
+    return python::make_tuple(a, dev, count);
 }
 
 using namespace boost::python;

src/graph/stats/graph_average.hh

@@ -27,6 +27,52 @@ namespace graph_tool
 using namespace std;
 using namespace boost;
 
+
+template <class Val1, class Val2>
+void operator+=(std::vector<Val1>& a, const std::vector<Val2>& b)
+{
+    a.resize(std::max(a.size(), b.size()));
+    for (size_t i = 0; i < std::min(a.size(), b.size()); ++i)
+        a[i] += b[i];
+}
+
+template <class Val1, class Val2>
+std::vector<Val1> operator*(const std::vector<Val1>& a, const std::vector<Val2>& b)
+{
+    std::vector<Val1> c(std::max(a.size(), b.size()));
+    for (size_t i = 0; i < std::min(a.size(), b.size()); ++i)
+        c[i] = a[i] * b[i];
+    return c;
+}
+
+template <class Val1, class Val2>
+std::vector<Val1> operator-(const std::vector<Val1>& a, const std::vector<Val2>& b)
+{
+    std::vector<Val1> c(std::max(a.size(), b.size()));
+    for (size_t i = 0; i < std::min(a.size(), b.size()); ++i)
+        c[i] = a[i] - b[i];
+    for (size_t i = a.size(); i < std::max(a.size(), b.size()); ++i)
+        c[i] = -b[i];
+    return c;
+}
+
+struct get_avg_type
+{
+    template <class Type>
+    struct apply
+    {
+        typedef typename mpl::if_<typename std::is_same<Type, python::object>::type,
+                                  python::object, long double>::type
+            type;
+    };
+
+    template <class Type>
+    struct apply<vector<Type>>
+    {
+        typedef vector<long double> type;
+    };
+};
+
 class VertexAverageTraverse
 {
 public:
@@ -35,9 +81,9 @@ public:
                     DegreeSelector& deg, ValueType& a, ValueType& aa,
                     size_t& count)
     {
-        ValueType x = deg(v, g);
+        const auto& x = deg(v, g);
         a += x;
-        aa += x*x;
+        aa += x * x;
         count++;
     }
 };
@@ -50,29 +96,48 @@ public:
                     EdgeProperty& eprop, ValueType& a, ValueType& aa,
                     size_t& count)
     {
-        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 (auto e : out_edges_range(v, g))
         {
-            ValueType x = eprop[*e];
+            const auto& x = eprop[e];
             a += x;
-            aa += x*x;
+            aa += x * x;
             count++;
         }
     }
 };
 
+// explicit special initialization to get around python::object
+template <class Val>
+void init_avg(Val& v)
+{
+    v = Val(0.);
+};
+
+template <class Val>
+void init_avg(std::vector<Val>&)
+{
+};
+
 // generalized functor to obtain average of different types of "degrees"
 template <class AverageTraverse>
 struct get_average
 {
-    get_average(long double& a, long double& dev)
-        : _a(a), _dev(dev) {}
+    get_average(boost::python::object& a, boost::python::object& dev,
+                size_t& count)
+        : _a(a), _dev(dev), _count(count) {}
 
     template <class Graph, class DegreeSelector>
     void operator()(Graph& g, DegreeSelector deg) const
     {
-        long double a = 0, aa = 0;
+        typedef typename DegreeSelector::value_type val_t;
+        dispatch(g, deg, typename std::is_pod<val_t>::type());
+    }
+
+    template <class Graph, class DegreeSelector>
+    void dispatch(Graph& g, DegreeSelector deg, std::true_type) const
+    {
+        typedef typename get_avg_type::apply<typename DegreeSelector::value_type>::type val_t;
+        val_t a = 0., aa = 0.;
         size_t count = 0;
 
         AverageTraverse traverse;
@@ -87,15 +152,34 @@ struct get_average
             traverse(g, v, deg, a, aa, count);
         }
 
-        _a = a/count;
-        _dev = sqrt((aa/count - _a*_a))/sqrt(count);
+        _a = boost::python::object(a);
+        _dev = boost::python::object(aa);
+        _count = count;
+    }
+
+    template <class Graph, class DegreeSelector>
+    void dispatch(Graph& g, DegreeSelector deg, std::false_type) const
+    {
+        typedef typename get_avg_type::apply<typename DegreeSelector::value_type>::type val_t;
+        val_t a, aa;
+        init_avg(a);
+        init_avg(aa);
+        size_t count = 0;
+
+        AverageTraverse traverse;
+        for (auto v : vertices_range(g))
+            traverse(g, v, deg, a, aa, count);
+
+        _a = boost::python::object(a);
+        _dev = boost::python::object(aa);
+        _count = count;
     }
 
-    long double& _a;
-    long double& _dev;
+    boost::python::object& _a;
+    boost::python::object& _dev;
+    size_t& _count;
 };
 
 } // graph_tool namespace
 
 #endif // GRAPH_AVERAGE_HH
-

src/graph_tool/stats/__init__.py

@@ -49,7 +49,7 @@ from __future__ import division, absolute_import, print_function
 from .. dl_import import dl_import
 dl_import("from . import libgraph_tool_stats")
 
-from .. import _degree, _prop, _get_rng, GraphView
+from .. import _degree, _prop, _get_rng, GraphView, PropertyMap
 from numpy import *
 import numpy
 import sys
@@ -234,9 +234,18 @@ def vertex_average(g, deg):
     (4.975, 0.0686758691244603)
     """
 
-    ret = libgraph_tool_stats.\
+    if isinstance(deg, PropertyMap) and "string" in deg.value_type():
+        raise ValueError("Cannot calculate average of property type: " + deg.value_type())
+    a, aa, count  = libgraph_tool_stats.\
           get_vertex_average(g._Graph__graph, _degree(g, deg))
-    return ret
+    try:
+        a = array(a.a)
+        aa = array(aa.a)
+    except AttributeError:
+        pass
+    a /= count
+    aa = sqrt((aa / count - a ** 2) / count)
+    return a, aa
 
 
 def edge_average(g, eprop):
@@ -287,9 +296,19 @@ def edge_average(g, eprop):
     (0.4989741369720412, 0.004101065927783255)
     """
 
-    ret = libgraph_tool_stats.\
+    if "string" in eprop.value_type():
+        raise ValueError("Cannot calculate average of property type: " + eprop.value_type())
+    g = GraphView(g, directed=True)
+    a, aa, count = libgraph_tool_stats.\
           get_edge_average(g._Graph__graph, _prop("e", g, eprop))
-    return ret
+    try:
+        a = array(a.a)
+        aa = array(aa.a)
+    except AttributeError:
+        pass
+    a /= count
+    aa = sqrt((aa / count - a ** 2) / count)
+    return a, aa
 
 
 def remove_labeled_edges(g, label):