graph_assortativity.hh 7.19 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2017 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
//
// 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 3
// 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, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_ASSORTATIVITY_HH
#define GRAPH_ASSORTATIVITY_HH

#include "shared_map.hh"
22
#include "graph_util.hh"
23
#include "hash_map_wrap.hh"
24

25
26
27
28
#if (BOOST_VERSION >= 106000)
# include <boost/math/special_functions/relative_difference.hpp>
#endif

29
30
31
32
33
namespace graph_tool
{
using namespace std;
using namespace boost;

Tiago Peixoto's avatar
Tiago Peixoto committed
34

35
36
37
38
39
40
// this will calculate the assortativity coefficient, based on the property
// pointed by 'deg'

struct get_assortativity_coefficient
{
    template <class Graph, class DegreeSelector>
41
    void operator()(const Graph& g, DegreeSelector deg, double& r,
42
43
                    double& r_err) const
    {
44
45
        size_t n_edges = 0;
        size_t e_kk = 0;
Tiago Peixoto's avatar
Tiago Peixoto committed
46
47

        typedef typename DegreeSelector::value_type val_t;
48
        typedef gt_hash_map<val_t, size_t> map_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
49
50
51
        map_t a, b;

        SharedMap<map_t> sa(a), sb(b);
52
53
        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
            firstprivate(sa, sb) reduction(+:e_kk, n_edges)
Tiago Peixoto's avatar
Tiago Peixoto committed
54
55
56
57
58
        parallel_vertex_loop_no_spawn
            (g,
             [&](auto v)
             {
                 val_t k1 = deg(v, g);
59
                 for (auto w : out_neighbors_range(v, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
60
61
62
                 {
                     val_t k2 = deg(w, g);
                     if (k1 == k2)
63
64
65
66
                         e_kk++;
                     sa[k1]++;
                     sb[k2]++;
                     n_edges++;
Tiago Peixoto's avatar
Tiago Peixoto committed
67
68
                 }
             });
69
70
71
72

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

Tiago Peixoto's avatar
Tiago Peixoto committed
73
        double t1 = double(e_kk) / n_edges, t2 = 0.0;
74

75
76
        for (auto& ai : a)
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
77
78
79
            auto bi = b.find(ai.first);
            if (bi != b.end())
                t2 += ai.second * bi->second;
80
        }
81
        t2 /= n_edges * n_edges;
82

83
84
85
86
87
88
89
90
91

#if (BOOST_VERSION >= 106000)
        if (boost::math::relative_difference(1., t2) > 1e-8)
#else
        if (abs(1.-t2) > 1e-8)
#endif
            r = (t1 - t2)/(1.0 - t2);
        else
            r = std::numeric_limits<double>::quiet_NaN();
92
93

        // "jackknife" variance
Tiago Peixoto's avatar
Tiago Peixoto committed
94
        double err = 0;
95
        size_t one = (graph_tool::is_directed(g)) ? 1 : 2;
96
97
        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
            reduction(+:err)
Tiago Peixoto's avatar
Tiago Peixoto committed
98
99
100
101
102
        parallel_vertex_loop_no_spawn
            (g,
             [&](auto v)
             {
                 val_t k1 = deg(v, g);
103
                 for (auto w : out_neighbors_range(v, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
104
105
                 {
                     val_t k2 = deg(w, g);
106
107
108
                     double tl2 = (t2 * (n_edges * n_edges)
                                   - one * b[k1] - one * a[k2]) /
                         ((n_edges - one) * (n_edges - one));
Tiago Peixoto's avatar
Tiago Peixoto committed
109
110
                     double tl1 = t1 * n_edges;
                     if (k1 == k2)
111
112
                         tl1 -= one;
                     tl1 /= n_edges - one;
Tiago Peixoto's avatar
Tiago Peixoto committed
113
                     double rl = (tl1 - tl2) / (1.0 - tl2);
114
                     err += (r - rl) * (r - rl);
Tiago Peixoto's avatar
Tiago Peixoto committed
115
116
                }
             });
117
        if (!graph_tool::is_directed(g))
118
            err /= 2;
119
120
121
122
123
124
125
126
#if (BOOST_VERSION >= 106000)
        if (boost::math::relative_difference(1., t2) > 1e-8)
#else
        if (abs(1.-t2) > 1e-8)
#endif
            r_err = sqrt(err);
        else
            r_err = std::numeric_limits<double>::quiet_NaN();
127
128
129
130
131
132
133
134
135
    }
};

// this will calculate the _scalar_ assortativity coefficient, based on the
// scalar property pointed by 'deg'

struct get_scalar_assortativity_coefficient
{
    template <class Graph, class DegreeSelector>
136
    void operator()(const Graph& g, DegreeSelector deg, double& r,
137
138
                    double& r_err) const
    {
139
        size_t n_edges = 0;
Tiago Peixoto's avatar
Tiago Peixoto committed
140
141
142
        double e_xy = 0;
        double a = 0, b = 0, da = 0, db = 0;

143
144
        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
            reduction(+:e_xy,n_edges,a,b,da,db)
Tiago Peixoto's avatar
Tiago Peixoto committed
145
146
147
148
        parallel_vertex_loop_no_spawn
            (g,
             [&](auto v)
             {
149
                 auto k1 = deg(v, g);
150
                 for (auto u : out_neighbors_range(v, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
151
152
                 {
                     auto k2 = deg(u, g);
153
154
155
156
157
158
                     a += k1;
                     da += k1 * k1;
                     b += k2;
                     db += k2 * k2;
                     e_xy += k1 * k2;
                     n_edges++;
Tiago Peixoto's avatar
Tiago Peixoto committed
159
160
                 }
             });
161
162

        double t1 = e_xy/n_edges;
163
164
        a /= n_edges;
        b /= n_edges;
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
        double stda;
        double stdb;

#if (BOOST_VERSION >= 106000)
        if (boost::math::relative_difference(da/n_edges, a*a) < 1e-8)
            stda = 0;
        else
            stda = sqrt(da/n_edges - a*a);
        if (boost::math::relative_difference(db/n_edges, b*b) < 1e-8)
            stdb = 0;
        else
            stdb = sqrt(db/n_edges - b*b);
#else
        if (sqrt(abs(da/n_edges - a*a)) < 1e-8)
            stda = 0;
        else
            stda = sqrt(da/n_edges - a*a);
        if (sqrt(abs(db/n_edges - b*b)) < 1e-8)
            stdb = 0;
        else
            stdb = sqrt(db/n_edges - b*b);
#endif
187

188
189
        if (stda*stdb > 0)
            r = (t1 - a*b)/(stda*stdb);
190
        else
191
            r = std::numeric_limits<double>::quiet_NaN();
192
193
194
195
196

        // "jackknife" variance
        r_err = 0.0;

        double err = 0.0;
197
        size_t one = (graph_tool::is_directed(g)) ? 1 : 2;
Tiago Peixoto's avatar
Tiago Peixoto committed
198
        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH) \
199
            reduction(+:err)
Tiago Peixoto's avatar
Tiago Peixoto committed
200
201
202
203
204
        parallel_vertex_loop_no_spawn
            (g,
             [&](auto v)
             {
                 double k1 = double(deg(v, g));
205
206
                 double al = (a * n_edges - k1) / (n_edges - one);
                 double dal = sqrt((da - k1 * k1) / (n_edges - one) - al * al);
Tiago Peixoto's avatar
Tiago Peixoto committed
207

208
                 for (auto u : out_neighbors_range(v, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
209
210
                 {
                     double k2 = deg(u, g);
211
212
213
                     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);
Tiago Peixoto's avatar
Tiago Peixoto committed
214
215
216
217
218
                     double rl;
                     if (dal * dbl > 0)
                         rl = (t1l - al * bl)/(dal * dbl);
                     else
                         rl = (t1l - al * bl);
219
                     err += (r - rl) * (r - rl);
Tiago Peixoto's avatar
Tiago Peixoto committed
220
221
                 }
             });
222
        if (!graph_tool::is_directed(g))
223
            err /= 2;
224
225
226
227
        if (stda*stdb > 0)
            r_err = sqrt(err);
        else
            r_err = std::numeric_limits<double>::quiet_NaN();
228
229
230
231
232
233
    }
};

} // graph_tool namespace

#endif //GRAPH_ASSORTATIVITY_HH