graph_assortativity.hh 6.81 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-2013 Tiago de Paula Peixoto <tiago@skewed.de>
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
//
// 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

21 22
#include "tr1_include.hh"
#include TR1_HEADER(unordered_map)
23

24
#include "shared_map.hh"
25
#include "graph_util.hh"
26 27 28 29 30 31 32 33 34 35 36 37

namespace graph_tool
{
using namespace std;
using namespace boost;

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

struct get_assortativity_coefficient
{
    template <class Graph, class DegreeSelector>
38
    void operator()(const Graph& g, DegreeSelector deg, double& r,
39 40
                    double& r_err) const
    {
41 42 43
        typedef typename mpl::if_<typename is_directed::apply<Graph>::type,
                                  size_t, double>::type count_t;

44
        count_t c = (is_directed::apply<Graph>::type::value) ? count_t(1) : count_t(0.5);
45 46 47 48
        count_t n_edges = 0;
        count_t e_kk = 0;
        tr1::unordered_map<double,count_t> a, b;
        SharedMap<tr1::unordered_map<double,count_t> > sa(a), sb(b);
49 50
        int i, N = num_vertices(g);
        #pragma omp parallel for default(shared) private(i) firstprivate(sa,sb)\
51
            schedule(static, 100) reduction(+:e_kk, n_edges)
52 53 54 55 56 57
        for (i = 0; i < N; ++i)
        {
            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
            if (v == graph_traits<Graph>::null_vertex())
                continue;

58
            double k1 = double(deg(v, g));
59 60 61
            typename graph_traits<Graph>::out_edge_iterator e, e_end;
            for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
            {
62
                double k2 = double(deg(target(*e, g), g));
63
                if (k1 == k2)
64 65 66 67
                    e_kk += c;
                sa[k1] += c;
                sb[k2] += c;
                n_edges += c;
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
            }
        }

        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)
            if (b.find(iter->second) != b.end())
                t2 += double(iter->second * b[iter->first]);
        t2 /= n_edges*n_edges;

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

        // "jackknife" variance
        double err = 0.0;
85
        #pragma omp parallel for default(shared) private(i) schedule(static, 100)\
86 87 88 89 90 91 92
            reduction(+:err)
        for (i = 0; i < N; ++i)
        {
            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
            if (v == graph_traits<Graph>::null_vertex())
                continue;

93
            double k1 = double(deg(v, g));
94 95 96
            typename graph_traits<Graph>::out_edge_iterator e, e_end;
            for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
            {
97
                double k2 = double(deg(target(*e,g),g));
98 99 100 101 102 103 104
                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);
105
                err += (r-rl)*(r-rl)*c;
106 107 108 109 110 111 112 113 114 115 116 117
            }
        }
        r_err = sqrt(err);
    }
};

// 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>
118
    void operator()(const Graph& g, DegreeSelector deg, double& r,
119 120
                    double& r_err) const
    {
121 122 123
        typedef typename mpl::if_<typename is_directed::apply<Graph>::type,
                                  size_t, double>::type count_t;

124
        count_t c = (is_directed::apply<Graph>::type::value) ? count_t(1) : count_t(0.5);
125
        count_t n_edges = 0;
126
        double e_xy = 0.0;
127
        double a = 0.0, b = 0.0, da = 0.0, db = 0.0;
128
        int i, N = num_vertices(g);
129
        #pragma omp parallel for default(shared) private(i) \
130
            schedule(static, 100) reduction(+:e_xy,n_edges,a,b,da,db)
131 132 133 134 135 136
        for (i = 0; i < N; ++i)
        {
            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
            if (v == graph_traits<Graph>::null_vertex())
                continue;

137
            double k1 = double(deg(v, g));
138 139 140
            typename graph_traits<Graph>::out_edge_iterator e, e_end;
            for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
            {
141
                double k2 = double(deg(target(*e,g),g));
142 143 144 145
                a += k1*c;
                da += k1*k1*c;
                b += k2*c;
                db += k2*k2*c;
146 147
                e_xy += k1*k2*c;
                n_edges += c;
148 149 150 151
            }
        }

        double t1 = e_xy/n_edges;
152 153 154 155
        a /= n_edges;
        b /= n_edges;
        double stda = sqrt(da/n_edges - a*a);
        double stdb = sqrt(db/n_edges - b*b);
156

157 158
        if (stda*stdb > 0)
            r = (t1 - a*b)/(stda*stdb);
159
        else
160
            r = (t1 - a*b);
161 162 163 164 165

        // "jackknife" variance
        r_err = 0.0;

        double err = 0.0;
166
        #pragma omp parallel for default(shared) private(i) schedule(static, 100)\
167 168 169 170 171 172 173
            reduction(+:err)
        for (i = 0; i < N; ++i)
        {
            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
            if (v == graph_traits<Graph>::null_vertex())
                continue;

174
            double k1 = double(deg(v, g));
175 176 177
            double al = (a*n_edges - k1)/(n_edges-1);
            double dal = sqrt((da - k1*k1)/(n_edges-1) - al*al);

178 179 180
            typename graph_traits<Graph>::out_edge_iterator e, e_end;
            for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
            {
181
                double k2 = double(deg(target(*e, g), g));
182 183
                double bl = (b*n_edges - k2)/(n_edges-1);
                double dbl = sqrt((db - k2*k2)/(n_edges-1) - bl*bl);
184 185 186
                double t1l = (e_xy - k1*k2)/(n_edges-1);
                double rl;
                if (dal*dbl > 0)
187
                    rl = (t1l - al*bl)/(dal*dbl);
188
                else
189
                    rl = (t1l - al*bl);
190
                err += (r-rl)*(r-rl)*c;
191 192 193 194 195 196 197 198 199
            }
        }
        r_err = sqrt(err);
    }
};

} // graph_tool namespace

#endif //GRAPH_ASSORTATIVITY_HH