* fix assortativity coefficient

* add jackknife variance
* add scalar assortativity coefficient


git-svn-id: https://svn.forked.de/graph-tool/trunk@11 d4600afd-f417-0410-95de-beed9576f240

src/graph-tool

@@ -105,6 +105,7 @@ correlations.add_option("--average-nearest-neighbours-correlation", action="call
 correlations.add_option("--edge-vertex-correlation-histogram", action="callback", callback=push_option, type="string", metavar="DEGREE1|EDGE-PROP|DEGREE2|FILE", help="get the source degree vs. edge scalar vs. target degree correlation histogram. Scalar properties are also accepted in place of DEGREE-TYPE")
 correlations.add_option("--vertex-scalar-correlation-histogram", action="callback", callback=push_option, type="string", metavar="DEGREE|VERTEX-PROP|FILE", help="get the degree vs. vertex scalar correlation histogram. Scalar properties are also accepted in place of DEGREE")
 correlations.add_option("--assortativity-coefficient", action="callback", callback=push_option, type="string", metavar="DEGREE|FILE", help="get the assortativity coefficient. Scalar properties are also accepted in place of DEGREE")
+correlations.add_option("--scalar-assortativity-coefficient", action="callback", callback=push_option, type="string", metavar="DEGREE|FILE", help="get the scalar assortativity coefficient. Scalar properties are also accepted in place of DEGREE")
 
 clustering =  parser.add_option_group("Clustering")
 clustering.add_option("--set-local-clustering-to-property", action="callback", callback=push_option, type="string", metavar="PROPERTY", help="set the local clustering coefficient to vertex property")
@@ -309,7 +310,19 @@ def parse_option(opt, just_file=False):
         deg = degree(values[0])
         if just_file:
             return values[1]
-        return (graph.GetAssortativityCoefficient(deg), values[1])
+        (r,err) = graph.GetAssortativityCoefficient(deg)
+        ret = "%f\t%f" % (r,err)
+        return (ret, values[1])
+    elif opt.name == "scalar-assortativity-coefficient":
+        values = parse_values(opt.value)
+        if len(values) != 2:
+            raise OptionError(opt.name, "invalid value '%s'" % opt.value)
+        deg = degree(values[0])
+        if just_file:
+            return values[1]
+        (r,err) = graph.GetScalarAssortativityCoefficient(deg)
+        ret = "%f\t%f" % (r,err)
+        return (ret, values[1])
     elif opt.name == "average-vertex-property" or opt.name == "average-edge-property":
         values = parse_values(opt.value)
         if len(values) != 2:

src/graph/graph.hh

@@ -78,7 +78,10 @@ public:
     hist3d_t   GetEdgeVertexCorrelationHistogram(deg_t deg1, std::string scalar, deg_t deg2) const;
     hist2d_t   GetVertexScalarCorrelationHistogram(deg_t deg, std::string scalar) const;
     avg_corr_t GetAverageNearestNeighboursCorrelation(deg_t origin_degree, deg_t neighbour_degree) const;
-    double     GetAssortativityCoefficient(deg_t deg) const;
+
+    // mixing
+    std::pair<double,double> GetAssortativityCoefficient(deg_t deg) const;
+    std::pair<double,double> GetScalarAssortativityCoefficient(deg_t deg) const;
 
     //clustering
     void   SetLocalClusteringToProperty(std::string property);

src/graph/graph_assortativity.cc

@@ -43,43 +43,62 @@ struct get_assortativity_coefficient
     get_assortativity_coefficient(DegreeSelector& deg): _deg(deg) {}
 
     template <class Graph>
-    void operator()(const Graph &g, double &a) const
+    void operator()(const Graph& g, double& r, double& r_err) const
     {
-	tr1::unordered_map<double,int> e_kk;
-	size_t n_edges = 0, n_vertices = 0;
-	tr1::unordered_map<double,int> n_k;
+	size_t n_edges = 0;
+	int e_kk = 0;
+	tr1::unordered_map<double,int> a, b;
 	
-	typename graph_traits<Graph>::vertex_iterator v, v_begin, v_end;
-	tie(v_begin,v_end) = vertices(g);
-	for(v = v_begin; v != v_end; ++v)
+	typename graph_traits<Graph>::edge_iterator e, e_begin, e_end;
+	tie(e_begin,e_end) = edges(g);
+	for (e = e_begin; e != e_end; ++e)
 	{
 	    double k1, k2;
-	    k1 = _deg(*v,g);
-	    n_k[k1]++;
-	    n_vertices++;
-	    
-	    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)
+	    k1 = _deg(source(*e,g),g);
+	    k2 = _deg(target(*e,g),g);
+	    a[k1]++;
+	    b[k2]++;
+	    if (k1 == k2)
+		e_kk++;
+	    n_edges++;
+	    if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
 	    {
-		k2 = _deg(target(*e,g),g);
+		a[k2]++;
+		a[k1]++;
 		if (k1 == k2)
-		    e_kk[k1]++;
+		    e_kk++;
 		n_edges++;
 	    }
 	}
 	
-	double t1=0.0, t2=0.0;
+	double t1=double(e_kk)/n_edges, t2=0.0;
 	
-	for (typeof(e_kk.begin()) iter = e_kk.begin(); iter != e_kk.end(); ++iter)
-	    t1 += double(iter->second)/n_edges;
+	for (typeof(a.begin()) iter = a.begin(); iter != a.end(); ++iter)
+	    if (b.find(iter->second) != b.end())
+		t2 += double(iter->second * b[iter->first]);
+	t2 /= n_edges*n_edges;
 	
-	double avg_k = GetHistogramMean(n_k);
-	
-	for (typeof(n_k.begin()) iter = n_k.begin(); iter != n_k.end(); ++iter)
-	    t2 += double(iter->second * iter->second * iter->first * iter->first)/(avg_k*avg_k*n_vertices*n_vertices);
-	
-	a = (t1 == 1.0)?1.0:(t1 - t2)/(1.0 - t2);
+	r = (t1 - t2)/(1.0 - t2);
+
+	// "jackknife" variance
+	r_err = 0.0;
+	for (e = e_begin; e != e_end; ++e)
+	{
+	    double k1, k2;
+	    k1 = _deg(source(*e,g),g);
+	    k2 = _deg(target(*e,g),g);
+	    int one = 1;
+	    if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
+		one = 2; // nice! :-)
+	    double tl2 = (t2*(n_edges*n_edges) - b[k1] - a[k2])/((n_edges-one)*(n_edges-one));
+	    double tl1 = t1*n_edges;
+	    if (k1==k2)
+		tl1 -= one;
+	    tl1 /= n_edges - one;
+	    double rl = (tl1 - tl2)/(1.0 - tl2);
+	    r_err += (r-rl)*(r-rl);
+	}
+	r_err = sqrt(r_err);
     }
     
     DegreeSelector& _deg;
@@ -87,8 +106,8 @@ struct get_assortativity_coefficient
 
 struct choose_assortativity_coefficient
 {
-    choose_assortativity_coefficient(const GraphInterface &g, GraphInterface::deg_t deg, double &a)
-	: _g(g), _a(a) 
+    choose_assortativity_coefficient(const GraphInterface& g, GraphInterface::deg_t deg, double& a, double& a_err)
+	: _g(g), _a(a), _a_err(a_err) 
     {
 	tie(_deg, _deg_name) = get_degree_type(deg);
     }
@@ -100,32 +119,157 @@ struct choose_assortativity_coefficient
 	if (mpl::at<degree_selector_index, DegreeSelector>::type::value == _deg)
 	{
 	    DegreeSelector deg(_deg_name, _g);
-	    check_filter(_g, bind<void>(get_assortativity_coefficient<DegreeSelector>(deg), _1, var(_a)), 
+	    check_filter(_g, bind<void>(get_assortativity_coefficient<DegreeSelector>(deg), _1, var(_a), var(_a_err)), 
 			 reverse_check(),directed_check());
 	}
     };
 
 
-    const GraphInterface &_g;
-    double &_a;
+    const GraphInterface& _g;
+    double& _a;
+    double& _a_err;
     GraphInterface::degree_t _deg;
     string _deg_name;
 };
 
 
-double 
+pair<double,double>
 GraphInterface::GetAssortativityCoefficient(GraphInterface::deg_t deg) const
 {
-    double a;
+    double a, a_err;
+    try
+    {
+	typedef mpl::vector<in_degreeS, out_degreeS, total_degreeS, scalarS> degrees;
+	mpl::for_each<degrees>(choose_assortativity_coefficient(*this, deg, a, a_err));
+    }
+    catch (dynamic_get_failure &e)
+    {
+	throw GraphException("error getting scalar property: " + string(e.what()));
+    }
+
+    return make_pair(a, a_err);
+}
+
+
+//==============================================================================
+// GetScalarAssortativityCoefficient(type)
+//==============================================================================
+
+template <class DegreeSelector>
+struct get_scalar_assortativity_coefficient
+{
+    get_scalar_assortativity_coefficient(DegreeSelector& deg): _deg(deg) {}
+
+    template <class Graph>
+    void operator()(const Graph& g, double& r, double& r_err) const
+    {
+	size_t n_edges = 0;
+	double e_xy = 0.0;
+	tr1::unordered_map<double,int> a, b;
+	
+	typename graph_traits<Graph>::edge_iterator e, e_begin, e_end;
+	tie(e_begin,e_end) = edges(g);
+	for (e = e_begin; e != e_end; ++e)
+	{
+	    double k1, k2;
+	    k1 = _deg(source(*e,g),g);
+	    k2 = _deg(target(*e,g),g);
+	    a[k1]++;
+	    b[k2]++;
+	    e_xy += k1*k2;
+	    n_edges++;
+	    if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
+	    {
+		a[k2]++;
+		b[k1]++;
+		e_xy += k1*k2;
+		n_edges++;
+	    }
+	}
+	
+	double t1 = e_xy/n_edges;
+	double avg_a = GetHistogramMean(a), avg_b = GetHistogramMean(b);
+	double sa = GetHistogramDeviation(a,avg_a), sb = GetHistogramDeviation(b,avg_b);
+	
+	if (sa*sb > 0)
+	    r = (t1 - avg_a*avg_b)/(sa*sb);
+	else
+	    r = (t1 - avg_a*avg_b);
+
+	// "jackknife" variance
+	r_err = 0.0;
+	double diff_a = 0.0, diff_b = 0.0;
+	for (typeof(a.begin()) iter = a.begin(); iter != a.end(); ++iter)
+	    diff_a += (iter->first - avg_a)*iter->second;
+	for (typeof(b.begin()) iter = b.begin(); iter != b.end(); ++iter)
+	    diff_b += (iter->first - avg_b)*iter->second;
+	for (e = e_begin; e != e_end; ++e)
+	{
+	    double k1, k2;
+	    k1 = _deg(source(*e,g),g);
+	    k2 = _deg(target(*e,g),g);
+	    int one = 1;
+	    if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
+		one = 2;
+	    double t1l = (e_xy - k1*k2)/(n_edges-one);
+	    double avg_al = (avg_a*n_edges - k1)/(n_edges-one), avg_bl = (avg_b*n_edges - k2)/(n_edges-one);
+	    double sal = sa - 2*diff_a*(avg_al-avg_a) + (avg_al-avg_a)*(avg_al-avg_a);
+	    double sbl = sb - 2*diff_b*(avg_bl-avg_b) + (avg_bl-avg_b)*(avg_bl-avg_b);
+	    double rl;
+	    if (sal*sbl > 0)
+		rl = (t1l - avg_al*avg_bl)/(sal*sbl);
+	    else
+		rl = (t1l - avg_al*avg_bl);
+	    r_err += (r-rl)*(r-rl);
+	}
+	r_err = sqrt(r_err);
+    }
+    
+    DegreeSelector& _deg;
+};
+
+struct choose_scalar_assortativity_coefficient
+{
+    choose_scalar_assortativity_coefficient(const GraphInterface& g, GraphInterface::deg_t deg, double& a, double& a_err)
+	: _g(g), _a(a), _a_err(a_err) 
+    {
+	tie(_deg, _deg_name) = get_degree_type(deg);
+    }
+    
+    
+    template <class DegreeSelector>
+    void operator()(DegreeSelector)
+    {	
+	if (mpl::at<degree_selector_index, DegreeSelector>::type::value == _deg)
+	{
+	    DegreeSelector deg(_deg_name, _g);
+	    check_filter(_g, bind<void>(get_scalar_assortativity_coefficient<DegreeSelector>(deg), _1, var(_a), var(_a_err)), 
+			 reverse_check(),directed_check());
+	}
+    };
+
+
+    const GraphInterface& _g;
+    double& _a;
+    double& _a_err;
+    GraphInterface::degree_t _deg;
+    string _deg_name;
+};
+
+
+pair<double,double>
+GraphInterface::GetScalarAssortativityCoefficient(GraphInterface::deg_t deg) const
+{
+    double a, a_err;
     try
     {
 	typedef mpl::vector<in_degreeS, out_degreeS, total_degreeS, scalarS> degrees;
-	mpl::for_each<degrees>(choose_assortativity_coefficient(*this, deg, a));
+	mpl::for_each<degrees>(choose_scalar_assortativity_coefficient(*this, deg, a, a_err));
     }
     catch (dynamic_get_failure &e)
     {
 	throw GraphException("error getting scalar property: " + string(e.what()));
     }
 
-    return a;
+    return make_pair(a, a_err);
 }

src/graph/graph_bind.cc

@@ -222,6 +222,7 @@ BOOST_PYTHON_MODULE(libgraph_tool)
 	.def("GetVertexScalarCorrelationHistogram", &GraphInterfaceWrap::GetVertexScalarCorrelationHistogram)
 	.def("GetAverageNearestNeighboursCorrelation", &GraphInterfaceWrap::GetAverageNearestNeighboursCorrelation)
 	.def("GetAssortativityCoefficient", &GraphInterfaceWrap::GetAssortativityCoefficient)
+	.def("GetScalarAssortativityCoefficient", &GraphInterfaceWrap::GetScalarAssortativityCoefficient)
 	.def("GetGlobalClustering", &GraphInterfaceWrap::GetGlobalClustering)
 	.def("SetLocalClusteringToProperty", &GraphInterfaceWrap::SetLocalClusteringToProperty)
 	.def("GetAverageDistance", &GraphInterfaceWrap::GetAverageDistance)

src/graph/histogram.hh

@@ -153,23 +153,38 @@ void ReadHistogram(Histogram &m, std::string input_file)
 
 //==============================================================================
 // GetHistogramMean()
-// Gets the mean value af a histogram
+// Gets the mean value of a histogram
 //==============================================================================
 template <class Histogram> 
 double GetHistogramMean (const Histogram &m)
 {
     int total = 0;
-    for (typeof(m.begin()) iter = m.begin(); iter != m.end(); iter++)
-        total += iter->second;
-  
-    if (total == 0) 
-        return 0.0;
-  
     double mean = 0;
     for (typeof(m.begin()) iter = m.begin(); iter != m.end(); iter++)
-        mean += double(iter->first * iter->second)/total;
+    {
+        mean += double(iter->first * iter->second);
+	total += iter->second;
+    }
     
-    return mean;
+    return (total > 0)?mean/total:0.0;
 }
 
+//==============================================================================
+// GetHistogramDeviation()
+// Gets the standard deviation of a histogram
+//==============================================================================
+template <class Histogram> 
+double GetHistogramDeviation (const Histogram &m, double avg)
+{  
+    double dev = 0.0;
+    int total = 0;
+    for (typeof(m.begin()) iter = m.begin(); iter != m.end(); iter++)
+    {
+        dev += double( (iter->first - avg) * (iter->first - avg) * iter->second);
+	total += iter->second;
+    }
+    return (total > 1)?sqrt(dev/(total-1)):0.0;
+}
+
+
 #endif //HISTOGRAM_HH