graph_assortativity.cc 10 KB
Newer Older
Tiago Peixoto's avatar
Tiago Peixoto committed
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
Tiago Peixoto's avatar
Tiago Peixoto committed
4
5
6
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
Tiago Peixoto's avatar
Tiago Peixoto committed
7
// as published by the Free Software Foundation; either version 3
Tiago Peixoto's avatar
Tiago Peixoto committed
8
9
10
11
12
13
14
15
// 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
16
17
// along with this program. If not, see <http://www.gnu.org/licenses/>.

Tiago Peixoto's avatar
Tiago Peixoto committed
18
19
20
21
22
23
24
25
26
27
28
29

#include <algorithm>
#include <tr1/unordered_set>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/random.hpp>

#include "graph.hh"
#include "histogram.hh"
#include "graph_filtering.hh"
#include "graph_selectors.hh"
#include "graph_properties.hh"
30
#include "shared_map.hh"
Tiago Peixoto's avatar
Tiago Peixoto committed
31
32
33
34
35
36

using namespace std;
using namespace boost;
using namespace boost::lambda;
using namespace graph_tool;

37
38
// this will calculate the assortativity coefficient, based on the property
// pointed by 'deg'
Tiago Peixoto's avatar
Tiago Peixoto committed
39

40
template <class DegreeSelector>
Tiago Peixoto's avatar
Tiago Peixoto committed
41
42
struct get_assortativity_coefficient
{
43
    get_assortativity_coefficient(DegreeSelector& deg): _deg(deg) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
44
45

    template <class Graph>
46
    void operator()(const Graph& g, double& r, double& r_err) const
Tiago Peixoto's avatar
Tiago Peixoto committed
47
    {
48
49
50
        size_t n_edges = 0;
        int e_kk = 0;
        tr1::unordered_map<double,int> a, b;
51
52
53
        SharedMap<tr1::unordered_map<double,int> > sa(a), sb(b);

        int i, N = num_vertices(g);
54
55
        #pragma omp parallel for default(shared) private(i) firstprivate(sa,sb)\
            schedule(dynamic) reduction(+:e_kk, n_edges)
56
        for (i = 0; i < N; ++i)
57
        {
58
59
60
61
62
63
64
            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)
65
            {
66
67
68
                double k2 = _deg(target(*e, g), g);
                sa[k1]++;
                sb[k2]++;
69
70
71
                if (k1 == k2)
                    e_kk++;
                n_edges++;
72
73
            }
        }
74
75
76
77

        sa.Gather();
        sb.Gather();

78
        double t1=double(e_kk)/n_edges, t2=0.0;
79

80
81
        for (typeof(a.begin()) iter = a.begin(); iter != a.end(); ++iter)
            if (b.find(iter->second) != b.end())
82
                t2 += double(iter->second * b[iter->first]);
83
        t2 /= n_edges*n_edges;
84

85
        r = (t1 - t2)/(1.0 - t2);
86

87
        // "jackknife" variance
88
        double err = 0.0;
89
90
        #pragma omp parallel for default(shared) private(i) schedule(dynamic)\
            reduction(+:err)
91
        for (i = 0; i < N; ++i)
92
        {
93
94
95
96
97
98
99
100
            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)
            {
101
102
103
                double k2 = _deg(target(*e,g),g);
                double tl2 = (t2*(n_edges*n_edges) - b[k1] - a[k2])/
                    ((n_edges-1)*(n_edges-1));
104
105
106
107
108
109
110
                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);
            }
111
        }
112
        r_err = sqrt(err);
Tiago Peixoto's avatar
Tiago Peixoto committed
113
    }
114

115
    DegreeSelector& _deg;
Tiago Peixoto's avatar
Tiago Peixoto committed
116
117
118
119
};

struct choose_assortativity_coefficient
{
120
121
122
123
    choose_assortativity_coefficient(const GraphInterface& g,
                                     GraphInterface::deg_t deg,
                                     double& a, double& a_err)
        : _g(g), _a(a), _a_err(a_err)
Tiago Peixoto's avatar
Tiago Peixoto committed
124
    {
125
        tie(_deg, _deg_name) = get_degree_type(deg);
Tiago Peixoto's avatar
Tiago Peixoto committed
126
    }
127
128


129
130
    template <class DegreeSelector>
    void operator()(DegreeSelector)
131
    {
132
133
        if (mpl::at<degree_selector_index, DegreeSelector>::type::value == _deg)
        {
134
            DegreeSelector deg(_deg_name, _g, true);
135
136
137
            check_filter(_g, bind<void>(get_assortativity_coefficient
                                            <DegreeSelector>(deg),
                                        _1, var(_a), var(_a_err)),
138
                         reverse_check(),directed_check());
139
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
140
141
    };

142
143
144
    const GraphInterface& _g;
    double& _a;
    double& _a_err;
145
146
    GraphInterface::degree_t _deg;
    string _deg_name;
Tiago Peixoto's avatar
Tiago Peixoto committed
147
148
149
};


150
pair<double,double>
151
GraphInterface::GetAssortativityCoefficient(GraphInterface::deg_t deg) const
Tiago Peixoto's avatar
Tiago Peixoto committed
152
{
153
154
155
    double a, a_err;
    try
    {
156
157
158
159
        typedef mpl::vector<in_degreeS, out_degreeS, total_degreeS, scalarS>
            degrees;
        mpl::for_each<degrees>(choose_assortativity_coefficient(*this, deg, a,
                                                                a_err));
160
161
162
    }
    catch (dynamic_get_failure &e)
    {
163
164
        throw GraphException("error getting scalar property: " +
                             string(e.what()));
165
166
167
168
169
    }

    return make_pair(a, a_err);
}

170
171
// this will calculate the _scalar_ assortativity coefficient, based on the
// scalar property pointed by 'deg'
172
173
174
175
176
177
178
179
180

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
    {
181
182
183
        size_t n_edges = 0;
        double e_xy = 0.0;
        tr1::unordered_map<double,int> a, b;
184
185
186
        SharedMap<tr1::unordered_map<double,int> > sa(a), sb(b);

        int i, N = num_vertices(g);
187
188
        #pragma omp parallel for default(shared) private(i) firstprivate(sa,sb)\
            schedule(dynamic) reduction(+:e_xy, n_edges)
189
        for (i = 0; i < N; ++i)
190
        {
191
192
193
194
195
196
197
            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)
198
            {
199
200
201
                double k2 = _deg(target(*e,g),g);
                sa[k1]++;
                sb[k2]++;
202
203
                e_xy += k1*k2;
                n_edges++;
204
205
            }
        }
206
207
208
209

        sa.Gather();
        sb.Gather();

210
211
        double t1 = e_xy/n_edges;
        double avg_a = GetHistogramMean(a), avg_b = GetHistogramMean(b);
212
213
214
        double da = GetHistogramDeviation(a,avg_a), db =
            GetHistogramDeviation(b,avg_b);

215
216
        if (da*db > 0)
            r = (t1 - avg_a*avg_b)/(da*db);
217
218
        else
            r = (t1 - avg_a*avg_b);
219

220
221
222
223
224
225
226
        // "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;
227

228
        double err = 0.0;
229
230
        #pragma omp parallel for default(shared) private(i) schedule(dynamic)\
            reduction(+:err)
231
        for (i = 0; i < N; ++i)
232
        {
233
234
235
236
237
238
239
240
            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)
            {
241
                double k2 = _deg(target(*e, g), g);
242
                double t1l = (e_xy - k1*k2)/(n_edges-1);
243
244
245
246
247
248
                double avg_al = (avg_a*n_edges - k1)/(n_edges-1);
                double 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);
249
250
251
252
253
254
255
                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);
            }
256
        }
257
        r_err = sqrt(err);
258
    }
259

260
261
262
263
264
    DegreeSelector& _deg;
};

struct choose_scalar_assortativity_coefficient
{
265
266
267
268
    choose_scalar_assortativity_coefficient(const GraphInterface& g,
                                            GraphInterface::deg_t deg,
                                            double& a, double& a_err)
        : _g(g), _a(a), _a_err(a_err)
269
    {
270
        tie(_deg, _deg_name) = get_degree_type(deg);
271
    }
272
273


274
275
    template <class DegreeSelector>
    void operator()(DegreeSelector)
276
    {
277
278
        if (mpl::at<degree_selector_index, DegreeSelector>::type::value == _deg)
        {
279
            DegreeSelector deg(_deg_name, _g, true);
280
281
282
            check_filter(_g, bind<void>(get_scalar_assortativity_coefficient
                                            <DegreeSelector>(deg),
                                        _1, var(_a), var(_a_err)),
283
                         reverse_check(),directed_check());
284
        }
285
286
287
288
289
290
291
292
293
294
295
    };

    const GraphInterface& _g;
    double& _a;
    double& _a_err;
    GraphInterface::degree_t _deg;
    string _deg_name;
};


pair<double,double>
296
297
GraphInterface::GetScalarAssortativityCoefficient(GraphInterface::deg_t deg)
    const
298
299
{
    double a, a_err;
Tiago Peixoto's avatar
Tiago Peixoto committed
300
301
    try
    {
302
303
304
305
        typedef mpl::vector<in_degreeS, out_degreeS, total_degreeS, scalarS>
            degrees;
        mpl::for_each<degrees>
            (choose_scalar_assortativity_coefficient(*this,deg, a, a_err));
Tiago Peixoto's avatar
Tiago Peixoto committed
306
307
308
    }
    catch (dynamic_get_failure &e)
    {
309
310
        throw GraphException("error getting scalar property: " +
                             string(e.what()));
Tiago Peixoto's avatar
Tiago Peixoto committed
311
312
    }

313
    return make_pair(a, a_err);
Tiago Peixoto's avatar
Tiago Peixoto committed
314
}