Commit 0c98262b authored by Tiago Peixoto's avatar Tiago Peixoto

* fixed bug in graph_assortativity.cc

* further code parellelization with openmp (vertex and edge histograms, correlations, and reciprocity)



git-svn-id: https://svn.forked.de/graph-tool/trunk@99 d4600afd-f417-0410-95de-beed9576f240
parent 0d0f4966
......@@ -100,14 +100,15 @@ AC_ARG_ENABLE([openmp], [AC_HELP_STRING([--enable-openmp],
if test $enableval = yes; then
[AC_MSG_RESULT(yes)]
[AC_DEFINE([USING_OPENMP], 1, [using openmp])]
USING_OPENMP=yes
[CXXFLAGS="$CXXFLAGS -fopenmp "]
USING_OPENMP=yes
[CXXFLAGS="$CXXFLAGS -fopenmp"]
[OPENMP_LDFLAGS=" -lgomp "]
fi
,
[AC_MSG_RESULT(no)]
[AC_DEFINE([USING_OPENMP], 0, [using openmp])]
USING_OPENMP=no
[CXXFLAGS="$CXXFLAGS -Wno-unknown-pragmas"]
[OPENMP_LDFLAGS=""]
)
AC_SUBST(OPENMP_LDFLAGS)
......
......@@ -44,6 +44,7 @@ libgraph_tool_la_SOURCES = \
graph_bind.cc\
graphml.cpp\
histogram.hh\
shared_map.hh\
read_graphviz_spirit.cpp
libgraph_tool_la_LIBADD = \
......
......@@ -39,6 +39,7 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
......@@ -232,10 +233,17 @@ struct get_vertex_histogram
template <class Graph, class Hist>
void operator()(const Graph& g, Hist& hist) const
{
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)
hist[_degree(*v,g)]++;
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[_degree(v,g)]++;
}
s_hist.Gather();
}
DegreeSelector& _degree;
};
......@@ -290,10 +298,25 @@ struct get_edge_histogram
template <class Graph, class Hist>
void operator()(const Graph& g, Hist& hist) const
{
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)
hist[_prop(*e,g)]++;
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;
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)
if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
s_hist[_prop(*e,g)] += 0.5;
else
s_hist[_prop(*e,g)]++;
}
s_hist.Gather();
}
scalarS& _prop;
};
......@@ -370,18 +393,23 @@ struct label_parallel_edges
template <class Graph, class EdgeIndexMap, class ParallelMap>
void operator()(const Graph& g, EdgeIndexMap edge_index, ParallelMap parallel) const
{
typename graph_traits<Graph>::vertex_iterator v, v_end;
for (tie(v, v_end) = vertices(g); v != v_end; ++v)
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;
tr1::unordered_set<typename graph_traits<Graph>::edge_descriptor,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)
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)
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);
......
......@@ -27,6 +27,7 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
......@@ -42,35 +43,39 @@ struct get_assortativity_coefficient
{
get_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;
int e_kk = 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)
SharedMap<tr1::unordered_map<double,int> > sa(a), sb(b);
int i, N = num_vertices(g);
#pragma omp parallel for default(shared) private(i) firstprivate(sa,sb) schedule(dynamic) reduction(+:e_kk, n_edges)
for (i = 0; i < N; ++i)
{
double k1, k2;
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)
typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
if (v == graph_traits<Graph>::null_vertex())
continue;
double k1 = _deg(v, g);
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
{
a[k2]++;
a[k1]++;
double k2 = _deg(target(*e, g), g);
sa[k1]++;
sb[k2]++;
if (k1 == k2)
e_kk++;
n_edges++;
}
}
sa.Gather();
sb.Gather();
double t1=double(e_kk)/n_edges, t2=0.0;
for (typeof(a.begin()) iter = a.begin(); iter != a.end(); ++iter)
......@@ -81,24 +86,29 @@ struct get_assortativity_coefficient
r = (t1 - t2)/(1.0 - t2);
// "jackknife" variance
r_err = 0.0;
for (e = e_begin; e != e_end; ++e)
double err = 0.0;
#pragma omp parallel for default(shared) private(i) schedule(dynamic) reduction(+:err)
for (i = 0; i < N; ++i)
{
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);
typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
if (v == graph_traits<Graph>::null_vertex())
continue;
double k1 = _deg(v, g);
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
{
double k2 = _deg(target(*e,g),g);
double tl2 = (t2*(n_edges*n_edges) - b[k1] - a[k2])/((n_edges-1)*(n_edges-1));
double tl1 = t1*n_edges;
if (k1 == k2)
tl1 -= 1;
tl1 /= n_edges - 1;
double rl = (tl1 - tl2)/(1.0 - tl2);
err += (r-rl)*(r-rl);
}
}
r_err = sqrt(r_err);
r_err = sqrt(err);
}
DegreeSelector& _deg;
......@@ -166,33 +176,38 @@ struct get_scalar_assortativity_coefficient
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)
SharedMap<tr1::unordered_map<double,int> > sa(a), sb(b);
int i, N = num_vertices(g);
#pragma omp parallel for default(shared) private(i) firstprivate(sa,sb) schedule(dynamic) reduction(+:e_xy, n_edges)
for (i = 0; i < N; ++i)
{
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)
typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
if (v == graph_traits<Graph>::null_vertex())
continue;
double k1 = _deg(v, g);
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
{
a[k2]++;
b[k1]++;
double k2 = _deg(target(*e,g),g);
sa[k1]++;
sb[k2]++;
e_xy += k1*k2;
n_edges++;
}
}
sa.Gather();
sb.Gather();
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);
double da = GetHistogramDeviation(a,avg_a), db = GetHistogramDeviation(b,avg_b);
if (sa*sb > 0)
r = (t1 - avg_a*avg_b)/(sa*sb);
if (da*db > 0)
r = (t1 - avg_a*avg_b)/(da*db);
else
r = (t1 - avg_a*avg_b);
......@@ -203,26 +218,33 @@ struct get_scalar_assortativity_coefficient
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 err = 0.0;
#pragma omp parallel for default(shared) private(i) schedule(dynamic) reduction(+:err)
for (i = 0; i < N; ++i)
{
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);
typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
if (v == graph_traits<Graph>::null_vertex())
continue;
double k1 = _deg(v, g);
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
{
double k2 = _deg(target(*e, g), g);
double t1l = (e_xy - k1*k2)/(n_edges-1);
double avg_al = (avg_a*n_edges - k1)/(n_edges-1), avg_bl = (avg_b*n_edges - k2)/(n_edges-1);
double dal = da - 2*diff_a*(avg_al-avg_a) + (avg_al-avg_a)*(avg_al-avg_a);
double dbl = db - 2*diff_b*(avg_bl-avg_b) + (avg_bl-avg_b)*(avg_bl-avg_b);
double rl;
if (dal*dbl > 0)
rl = (t1l - avg_al*avg_bl)/(dal*dbl);
else
rl = (t1l - avg_al*avg_bl);
err += (r-rl)*(r-rl);
}
}
r_err = sqrt(r_err);
r_err = sqrt(err);
}
DegreeSelector& _deg;
......
......@@ -27,6 +27,7 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
......@@ -44,22 +45,28 @@ struct get_correlation_histogram
get_correlation_histogram(DegreeSelector1& deg1, DegreeSelector2& deg2): _deg1(deg1), _deg2(deg2) {}
template <class Graph, class WeightMap, class Hist>
void operator()(Graph &g, WeightMap weight, Hist &hist) const
void operator()(Graph& g, WeightMap weight, Hist& hist) const
{
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)
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;
typename Hist::key_type key;
key.first = _deg1(source(*e,g),g);
key.second = _deg2(target(*e,g),g);
hist[key] += typename Hist::value_type::second_type(get(weight, *e));
if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
key.first = _deg1(v, g);
typename graph_traits<Graph>::out_edge_iterator e, e_end;
for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
{
swap(key.first, key.second);
hist[key] += get(weight, *e);
key.second = _deg2(target(*e,g),g);
s_hist[key] += typename Hist::value_type::second_type(get(weight, *e));
}
}
s_hist.Gather();
}
DegreeSelector1& _deg1;
DegreeSelector2& _deg2;
......@@ -155,5 +162,3 @@ GraphInterface::GetVertexCorrelationHistogram(GraphInterface::deg_t deg1, GraphI
return hist;
}
......@@ -27,6 +27,7 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
......@@ -46,13 +47,23 @@ struct get_combined_degree_histogram
template <class Graph, class Hist>
void operator()(Graph &g, Hist &hist) const
{
SharedMap<Hist> s_hist(hist);
typedef typename Hist::key_type::first_type first_type;
typedef typename Hist::key_type::second_type second_type;
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)
hist[make_pair(first_type(_deg1(*v,g)), second_type(_deg2(*v,g)))]++;
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[make_pair(first_type(_deg1(v,g)), second_type(_deg2(v,g)))]++;
}
s_hist.Gather();
}
DegreeSelector1& _deg1;
DegreeSelector2& _deg2;
......
......@@ -27,12 +27,19 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
using namespace boost::lambda;
using namespace graph_tool;
void operator+=(pair<double, double>&a, const pair<double, double>&b)
{
a.first += b.first;
a.second += b.second;
}
//==============================================================================
// average_nearest_neighbours_correlation
// return generalized average nearest neighbours correlation
......@@ -48,23 +55,33 @@ struct get_average_nearest_neighbours_correlation
void operator()(const Graph& g, WeightMap weight, AvgDeg& avg_deg) const
{
tr1::unordered_map<double,double> count;
SharedMap<tr1::unordered_map<double,double> > s_count(count);
SharedMap<AvgDeg> s_avg_deg(avg_deg);
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)
int i, N = num_vertices(g);
#pragma omp parallel for default(shared) private(i) firstprivate(s_count, s_avg_deg) 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;
typename AvgDeg::key_type orig_deg = _origin_degree(v,g);
typename graph_traits<Graph>::out_edge_iterator e, e_begin, e_end;
tie(e_begin,e_end) = out_edges(*v,g);
tie(e_begin,e_end) = out_edges(v,g);
for(e = e_begin; e != e_end; ++e)
{
{
typename AvgDeg::value_type::second_type::first_type deg = _neighbours_degree(target(*e,g),g);
typename AvgDeg::key_type orig_deg = _origin_degree(*v,g);
avg_deg[orig_deg].first += deg*get(weight, *e);
avg_deg[orig_deg].second += deg*deg;
count[orig_deg] += get(weight,*e);
s_avg_deg[orig_deg].first += deg*get(weight, *e);
s_avg_deg[orig_deg].second += deg*deg;
s_count[orig_deg] += get(weight,*e);
}
}
s_count.Gather();
s_avg_deg.Gather();
for (typeof(avg_deg.begin()) iter = avg_deg.begin(); iter != avg_deg.end(); ++iter)
{
double N = count[iter->first];
......
......@@ -27,6 +27,7 @@
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
#include "shared_map.hh"
using namespace std;
using namespace boost;
......@@ -47,21 +48,29 @@ struct get_edge_correlation_histogram
template <class Graph, class Hist>
void operator()(Graph &g, Hist &hist) const
{
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)
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;
typename Hist::key_type key;
get<0>(key) = _deg1(source(*e,g),g);
get<1>(key) = _edge_scalar(*e, g);
get<2>(key) = _deg2(target(*e,g),g);
hist[key]++;
if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
get<0>(key) = _deg1(v,g);
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)
{
swap(get<0>(key), get<2>(key));
hist[key]++;
get<1>(key) = _edge_scalar(*e, g);
get<2>(key) = _deg2(target(*e,g),g);
s_hist[key]++;
}
}
s_hist.Gather();
}
scalarS& _edge_scalar;
DegreeSelector1& _deg1;
......
......@@ -38,26 +38,42 @@ struct get_reciprocity
size_t L = 0;
double Lbd = 0.0;
typename graph_traits<Graph>::edge_iterator e,e_end;
for (tie(e,e_end) = edges(g); e != e_end; ++e)
int i, NV = num_vertices(g);
#pragma omp parallel for default(shared) private(i) reduce(+:L,Ldb) schedule(dynamic)
for (i = 0; i < NV; ++i)
{
typename graph_traits<Graph>::vertex_descriptor s,t;
s = source(*e,g);
t = target(*e,g);
typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
if (v == graph_traits<Graph>::null_vertex())
continue;
size_t o = 0;
typename graph_traits<Graph>::adjacency_iterator a, a_end;
for (tie(a,a_end) = adjacent_vertices(s, g); a != a_end; ++a)
if (*a == t)
o++;
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 s,t;
s = v;
t = target(*e,g);
size_t i = 0;
for (tie(a, a_end) = adjacent_vertices(t, g); a != a_end; ++a)
if (*a == s)
i++;
size_t o = 0;
typename graph_traits<Graph>::adjacency_iterator a, a_end;
for (tie(a,a_end) = adjacent_vertices(s, g); a != a_end; ++a)
if (*a == t)
o++;
Lbd += min(i/double(o),1.0);
L++;
size_t j = 0;
for (tie(a, a_end) = adjacent_vertices(t, g); a != a_end; ++a)
if (*a == s)
j++;
Lbd += min(j/double(o),1.0);
L++;
}
}
if(is_convertible<typename graph_traits<Graph>::directed_category, undirected_tag>::value)
{
L /= 2;
Lbd /= 2;
}
size_t N = HardNumVertices()(g);
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006 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 2
// 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, write to the Free Software