Commit 1495bc72 authored by Tiago Peixoto's avatar Tiago Peixoto

Implement weighted clustering and assortativity

parent 81232cbf
......@@ -31,22 +31,45 @@ using namespace std;
using namespace boost;
using namespace graph_tool;
boost::python::tuple global_clustering(GraphInterface& g)
boost::python::tuple global_clustering(GraphInterface& g, boost::any weight)
{
typedef UnityPropertyMap<size_t,GraphInterface::edge_t> weight_map_t;
typedef boost::mpl::push_back<edge_scalar_properties, weight_map_t>::type
weight_props_t;
if (!weight.empty() && !belongs<edge_scalar_properties>()(weight))
throw ValueException("weight edge property must have a scalar value type");
if(weight.empty())
weight = weight_map_t();
double c, c_err;
run_action<graph_tool::detail::never_directed>()
(g, std::bind(get_global_clustering(), std::placeholders::_1,
std::ref(c), std::ref(c_err)))();
std::placeholders::_2, std::ref(c), std::ref(c_err)),
weight_props_t())(weight);
return boost::python::make_tuple(c, c_err);
}
void local_clustering(GraphInterface& g, boost::any prop)
void local_clustering(GraphInterface& g, boost::any prop, boost::any weight)
{
typedef UnityPropertyMap<size_t,GraphInterface::edge_t> weight_map_t;
typedef boost::mpl::push_back<edge_scalar_properties, weight_map_t>::type
weight_props_t;
if (!weight.empty() && !belongs<edge_scalar_properties>()(weight))
throw ValueException("weight edge property must have a scalar value type");
if(weight.empty())
weight = weight_map_t();
run_action<>()
(g, std::bind(set_clustering_to_property(),
std::placeholders::_1,
std::placeholders::_2),
writable_vertex_scalar_properties())(prop);
std::placeholders::_2,
std::placeholders::_3),
weight_props_t(),
writable_vertex_scalar_properties())(weight, prop);
}
using namespace boost::python;
......
......@@ -44,52 +44,58 @@ using namespace boost;
using namespace std;
// calculates the number of triangles to which v belongs
template <class Graph, class VProp>
pair<int,int>
get_triangles(typename graph_traits<Graph>::vertex_descriptor v, VProp& mark,
const Graph& g)
template <class Graph, class EWeight, class VProp>
auto get_triangles(typename graph_traits<Graph>::vertex_descriptor v,
EWeight& eweight, VProp& mark, const Graph& g)
{
size_t triangles = 0;
typedef typename property_traits<EWeight>::value_type val_t;
val_t triangles = 0, w1 = 0, w2 = 0, k = 0;
for (auto n : adjacent_vertices_range(v, g))
for (auto e : out_edges_range(v, g))
{
auto n = target(e, g);
if (n == v)
continue;
mark[n] = true;
auto w = eweight[e];
mark[n] = w;
k += w;
}
for (auto n : adjacent_vertices_range(v, g))
for (auto e : out_edges_range(v, g))
{
auto n = target(e, g);
if (n == v)
continue;
for (auto n2 : adjacent_vertices_range(n, g))
w1 = eweight[e];
for (auto e2 : out_edges_range(n, g))
{
auto n2 = target(e2, g);
if (n2 == n)
continue;
if (mark[n2])
++triangles;
w2 = eweight[e2];
triangles += mark[n2] * w1 * w2;
}
}
for (auto n : adjacent_vertices_range(v, g))
mark[n] = false;
mark[n] = 0;
size_t k = out_degree(v, g);
if (graph_tool::is_directed(g))
return make_pair(triangles, (k * (k - 1)));
return make_pair(val_t(triangles), val_t((k * (k - 1))));
else
return make_pair(triangles / 2, (k * (k - 1)) / 2);
return make_pair(val_t(triangles / 2), val_t((k * (k - 1)) / 2));
}
// retrieves the global clustering coefficient
struct get_global_clustering
{
template <class Graph>
void operator()(const Graph& g, double& c, double& c_err) const
template <class Graph, class EWeight>
void operator()(const Graph& g, EWeight eweight, double& c, double& c_err) const
{
size_t triangles = 0, n = 0;
vector<bool> mask(num_vertices(g), false);
typedef typename property_traits<EWeight>::value_type val_t;
val_t triangles = 0, n = 0;
vector<val_t> mask(num_vertices(g), 0);
#pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
firstprivate(mask) reduction(+:triangles, n)
......@@ -97,7 +103,7 @@ struct get_global_clustering
(g,
[&](auto v)
{
auto temp = get_triangles(v, mask, g);
auto temp = get_triangles(v, eweight, mask, g);
triangles += temp.first;
n += temp.second;
});
......@@ -112,7 +118,7 @@ struct get_global_clustering
(g,
[&](auto v)
{
auto temp = get_triangles(v, mask, g);
auto temp = get_triangles(v, eweight, mask, g);
double cl = double(triangles - temp.first) /
(n - temp.second);
cerr += power(c - cl, 2);
......@@ -124,11 +130,11 @@ struct get_global_clustering
// sets the local clustering coefficient to a property
struct set_clustering_to_property
{
template <class Graph, class ClustMap>
void operator()(const Graph& g, ClustMap clust_map) const
template <class Graph, class EWeight, class ClustMap>
void operator()(const Graph& g, EWeight eweight, ClustMap clust_map) const
{
typedef typename property_traits<ClustMap>::value_type c_type;
vector<bool> mask(num_vertices(g), false);
typedef typename property_traits<EWeight>::value_type val_t;
vector<val_t> mask(num_vertices(g), false);
#pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
firstprivate(mask)
......@@ -136,11 +142,11 @@ struct set_clustering_to_property
(g,
[&](auto v)
{
auto triangles = get_triangles(v, mask, g);
auto triangles = get_triangles(v, eweight, mask, g);
double clustering = (triangles.second > 0) ?
double(triangles.first)/triangles.second :
0.0;
clust_map[v] = c_type(clustering);
clust_map[v] = clustering;
});
}
......
......@@ -26,28 +26,50 @@ using namespace std;
using namespace graph_tool;
pair<double,double>
assortativity_coefficient(GraphInterface& gi,
GraphInterface::deg_t deg)
assortativity_coefficient(GraphInterface& gi, GraphInterface::deg_t deg,
boost::any weight)
{
typedef UnityPropertyMap<size_t,GraphInterface::edge_t> weight_map_t;
typedef boost::mpl::push_back<edge_scalar_properties, weight_map_t>::type
weight_props_t;
if (!weight.empty() && !belongs<edge_scalar_properties>()(weight))
throw ValueException("weight edge property must have a scalar value type");
if(weight.empty())
weight = weight_map_t();
double a = 0, a_err = 0;
run_action<>()(gi,std::bind(get_assortativity_coefficient(),
std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3,
std::ref(a), std::ref(a_err)),
scalar_selectors())
(degree_selector(deg));
scalar_selectors(), weight_props_t())
(degree_selector(deg), weight);
return make_pair(a, a_err);
}
pair<double,double>
scalar_assortativity_coefficient(GraphInterface& gi,
GraphInterface::deg_t deg)
scalar_assortativity_coefficient(GraphInterface& gi, GraphInterface::deg_t deg,
boost::any weight)
{
typedef UnityPropertyMap<size_t,GraphInterface::edge_t> weight_map_t;
typedef boost::mpl::push_back<edge_scalar_properties, weight_map_t>::type
weight_props_t;
if (!weight.empty() && !belongs<edge_scalar_properties>()(weight))
throw ValueException("weight edge property must have a scalar value type");
if(weight.empty())
weight = weight_map_t();
double a = 0, a_err = 0;
run_action<>()(gi, std::bind(get_scalar_assortativity_coefficient(),
std::placeholders::_1, std::placeholders::_2,
std::placeholders::_3,
std::ref(a), std::ref(a_err)),
scalar_selectors())
(degree_selector(deg));
scalar_selectors(), weight_props_t())
(degree_selector(deg), weight);
return make_pair(a, a_err);
}
......
......@@ -37,12 +37,13 @@ using namespace boost;
struct get_assortativity_coefficient
{
template <class Graph, class DegreeSelector>
void operator()(const Graph& g, DegreeSelector deg, double& r,
double& r_err) const
template <class Graph, class DegreeSelector, class Eweight>
void operator()(const Graph& g, DegreeSelector deg, Eweight eweight,
double& r, double& r_err) const
{
size_t n_edges = 0;
size_t e_kk = 0;
typedef typename property_traits<Eweight>::value_type wval_t;
wval_t n_edges = 0;
wval_t e_kk = 0;
typedef typename DegreeSelector::value_type val_t;
typedef gt_hash_map<val_t, size_t> map_t;
......@@ -56,14 +57,16 @@ struct get_assortativity_coefficient
[&](auto v)
{
val_t k1 = deg(v, g);
for (auto w : out_neighbors_range(v, g))
for (auto e : out_edges_range(v, g))
{
val_t k2 = deg(w, g);
auto u = target(e, g);
auto w = eweight[e];
val_t k2 = deg(u, g);
if (k1 == k2)
e_kk++;
sa[k1]++;
sb[k2]++;
n_edges++;
e_kk += w;
sa[k1] += w;
sb[k2] += w;
n_edges += w;
}
});
......@@ -100,16 +103,18 @@ struct get_assortativity_coefficient
[&](auto v)
{
val_t k1 = deg(v, g);
for (auto w : out_neighbors_range(v, g))
for (auto e : out_edges_range(v, g))
{
val_t k2 = deg(w, g);
auto u = target(e, g);
auto w = eweight[e];
val_t k2 = deg(u, g);
double tl2 = (t2 * (n_edges * n_edges)
- one * b[k1] - one * a[k2]) /
((n_edges - one) * (n_edges - one));
- w * one * b[k1] - w * one * a[k2]) /
((n_edges - w * one) * (n_edges - w * one));
double tl1 = t1 * n_edges;
if (k1 == k2)
tl1 -= one;
tl1 /= n_edges - one;
tl1 -= one * w;
tl1 /= n_edges - one * w;
double rl = (tl1 - tl2) / (1.0 - tl2);
err += (r - rl) * (r - rl);
}
......@@ -132,11 +137,12 @@ struct get_assortativity_coefficient
struct get_scalar_assortativity_coefficient
{
template <class Graph, class DegreeSelector>
void operator()(const Graph& g, DegreeSelector deg, double& r,
double& r_err) const
template <class Graph, class DegreeSelector, class Eweight>
void operator()(const Graph& g, DegreeSelector deg, Eweight eweight,
double& r, double& r_err) const
{
size_t n_edges = 0;
typedef typename property_traits<Eweight>::value_type val_t;
val_t n_edges = 0;
double e_xy = 0;
double a = 0, b = 0, da = 0, db = 0;
......@@ -147,15 +153,17 @@ struct get_scalar_assortativity_coefficient
[&](auto v)
{
auto k1 = deg(v, g);
for (auto u : out_neighbors_range(v, g))
for (auto e : out_edges_range(v, g))
{
auto u = target(e, g);
auto w = eweight[e];
auto k2 = deg(u, g);
a += k1;
da += k1 * k1;
b += k2;
db += k2 * k2;
e_xy += k1 * k2;
n_edges++;
a += k1 * w;
da += k1 * k1 * w;
b += k2 * w;
db += k2 * k2 * w;
e_xy += k1 * k2 * w;
n_edges += w;
}
});
......@@ -205,12 +213,14 @@ struct get_scalar_assortativity_coefficient
double al = (a * n_edges - k1) / (n_edges - one);
double dal = sqrt((da - k1 * k1) / (n_edges - one) - al * al);
for (auto u : out_neighbors_range(v, g))
for (auto e : out_edges_range(v, g))
{
auto u = target(e, g);
auto w = eweight[e];
double k2 = deg(u, g);
double bl = (b * n_edges - k2 * one) / (n_edges - one);
double dbl = sqrt((db - k2 * k2 * one) / (n_edges - one) - bl * bl);
double t1l = (e_xy - k1 * k2 * one)/(n_edges - one);
double bl = (b * n_edges - k2 * one * w) / (n_edges - one * w);
double dbl = sqrt((db - k2 * k2 * one * w) / (n_edges - one * w) - bl * bl);
double t1l = (e_xy - k1 * k2 * one * w)/(n_edges - one * w);
double rl;
if (dal * dbl > 0)
rl = (t1l - al * bl)/(dal * dbl);
......
......@@ -60,7 +60,7 @@ __all__ = ["local_clustering", "global_clustering", "extended_clustering",
"motifs", "motif_significance"]
def local_clustering(g, prop=None, undirected=True):
def local_clustering(g, weight=None, prop=None, undirected=True):
r"""
Return the local clustering coefficients for all vertices.
......@@ -68,6 +68,8 @@ def local_clustering(g, prop=None, undirected=True):
----------
g : :class:`~graph_tool.Graph`
Graph to be used.
weight : :class:`~graph_tool.PropertyMap`, optional (default: None)
Edge weights. If omitted, a constant value of 1 will be used.
prop : :class:`~graph_tool.PropertyMap` or string, optional
Vertex property map where results will be stored. If specified, this
parameter will also be the return value.
......@@ -105,8 +107,8 @@ def local_clustering(g, prop=None, undirected=True):
.. math::
c'_i = 2c_i.
The implemented algorithm runs in :math:`O(|V|\left< k\right>^2)` time,
where :math:`\left< k\right>` is the average out-degree.
The implemented algorithm runs in :math:`O(|V|\left<k^2\right>)` time,
where :math:`\left<k^2\right>` is second moment of the degree distribution.
If enabled during compilation, this algorithm runs in parallel.
......@@ -128,18 +130,20 @@ def local_clustering(g, prop=None, undirected=True):
prop = g.new_vertex_property("double")
if g.is_directed() and undirected:
g = GraphView(g, directed=False, skip_properties=True)
_gt.local_clustering(g._Graph__graph, _prop("v", g, prop))
_gt.local_clustering(g._Graph__graph, _prop("v", g, prop),
_prop("e", g, weight))
return prop
def global_clustering(g):
r"""
Return the global clustering coefficient.
def global_clustering(g, weight=None):
r"""Return the global clustering coefficient.
Parameters
----------
g : :class:`~graph_tool.Graph`
Graph to be used.
weight : :class:`~graph_tool.PropertyMap`, optional (default: None)
Edge weights. If omitted, a constant value of 1 will be used.
Returns
-------
......@@ -161,8 +165,17 @@ def global_clustering(g):
c = 3 \times \frac{\text{number of triangles}}
{\text{number of connected triples}}
The implemented algorithm runs in :math:`O(|V|\left< k\right>^2)` time,
where :math:`\left< k\right>` is the average (total) degree.
If weights are given, the following definition is used
.. math::
c = \frac{\operatorname{Tr}{{\boldsymbol A}^3}}{\sum_{i\ne j}[{\boldsymbol A}^2]_{ij}},
where :math:`\boldsymbol A` is the weighted adjacency matrix, and it is
assumed that the weights are normalized, i.e. :math:`A_{ij} \le 1`.
The implemented algorithm runs in time :math:`O(|V|\left<k^2\right>)`,
where :math:`\left< k^2\right>` is the second moment of the degree
distribution.
If enabled during compilation, this algorithm runs in parallel.
......@@ -177,11 +190,12 @@ def global_clustering(g):
.. [newman-structure-2003] M. E. J. Newman, "The structure and function of
complex networks", SIAM Review, vol. 45, pp. 167-256, 2003,
:doi:`10.1137/S003614450342480`
"""
if g.is_directed():
g = GraphView(g, directed=False, skip_properties=True)
c = _gt.global_clustering(g._Graph__graph)
c = _gt.global_clustering(g._Graph__graph, _prop("e", g, weight))
return c
......
......@@ -52,17 +52,19 @@ __all__ = ["assortativity", "scalar_assortativity", "corr_hist",
"avg_combined_corr"]
def assortativity(g, deg):
r"""
Obtain the assortativity coefficient for the given graph.
def assortativity(g, deg, eweight=None):
r"""Obtain the assortativity coefficient for the given graph.
Parameters
----------
g : :class:`~graph_tool.Graph`
Graph to be used.
deg : string or :class:`~graph_tool.PropertyMap`
degree type ("in", "out" or "total") or vertex property map, which
Degree type ("in", "out" or "total") or vertex property map, which
specifies the vertex types.
eweight : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
If given, this will specify the edge weights, otherwise a constant value
of one will be used.
Returns
-------
......@@ -114,10 +116,11 @@ def assortativity(g, deg):
"""
return libgraph_tool_correlations.\
assortativity_coefficient(g._Graph__graph, _degree(g, deg))
assortativity_coefficient(g._Graph__graph, _degree(g, deg),
_prop("e", g, eweight))
def scalar_assortativity(g, deg):
def scalar_assortativity(g, deg, eweight=None):
r"""
Obtain the scalar assortativity coefficient for the given graph.
......@@ -126,8 +129,11 @@ def scalar_assortativity(g, deg):
g : :class:`~graph_tool.Graph`
Graph to be used.
deg : string or :class:`~graph_tool.PropertyMap`
degree type ("in", "out" or "total") or vertex property map, which
specifies the vertex types.
Degree type ("in", "out" or "total") or vertex property map, which
specifies the vertex scalar values.
eweight : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
If given, this will specify the edge weights, otherwise a constant value
of one will be used.
Returns
-------
......@@ -177,8 +183,8 @@ def scalar_assortativity(g, deg):
.. _jackknife method: http://en.wikipedia.org/wiki/Resampling_%28statistics%29#Jackknife
"""
return libgraph_tool_correlations.\
scalar_assortativity_coefficient(g._Graph__graph,
_degree(g, deg))
scalar_assortativity_coefficient(g._Graph__graph, _degree(g, deg),
_prop("e", g, eweight))
def corr_hist(g, deg_source, deg_target, bins=[[0, 1], [0, 1]], weight=None,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment