graph_correlations.hh 7.12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2007  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 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_CORRELATIONS_HH
#define GRAPH_CORRELATIONS_HH

#include <algorithm>
22
23
24
25
26
27
#include <boost/numeric/conversion/bounds.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/python/object.hpp>
#include <boost/python/list.hpp>
#include <boost/python/extract.hpp>
#include "numpy_bind.hh"
28
29
30
31
32
33
34
35
#include "shared_map.hh"

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

36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
// get degrees pairs from source and of neighbours
class GetNeighboursPairs
{
public:

    template <class Graph, class Deg1, class Deg2, class Hist, class WeightMap>
    void operator()(typename graph_traits<Graph>::vertex_descriptor v,
                    Deg1& deg1, Deg2& deg2, Graph& g, WeightMap& weight,
                    Hist& hist)
    {
        typename Hist::point_t k;
        k[0] = 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)
        {
            k[1] = deg2(target(*e,g),g);
            hist.PutValue(k, get(weight, *e));
        }
    }
};

// get degrees pairs from one single vertex
class GetCombinedPair
{
public:

    template <class Graph, class Deg1, class Deg2, class Hist, class Dummy>
    void operator()(typename graph_traits<Graph>::vertex_descriptor v,
                    Deg1& deg1, Deg2& deg2, Graph& g, const Dummy&,
                    Hist& hist)
    {
        typename Hist::point_t k;
        k[0] = deg1(v, g);
        k[1] = deg2(v, g);
        hist.PutValue(k);
    }
};


namespace detail
{
struct select_larger_type
{
    template <class Type1, class Type2>
    struct apply
    {
        typedef typename mpl::if_<
            typename mpl::greater<typename mpl::sizeof_<Type1>::type,
                                  typename mpl::sizeof_<Type2>::type>::type,
            Type1,
            Type2>::type type;
    };
};

struct select_float_and_larger
{
    template <class Type1, class Type2>
    struct apply
    {
        typedef typename mpl::if_<
            mpl::equal_to<is_floating_point<Type1>,
                          is_floating_point<Type2> >,
            typename select_larger_type::apply<Type1,Type2>::type,
            typename mpl::if_<is_floating_point<Type1>,
                              Type1,
                              Type2>::type>::type type;
    };
};

}
106

107
108
// retrieves the generalized vertex-vertex correlation histogram
template <class GetDegreePair>
109
110
struct get_correlation_histogram
{
111
112
113
114
    get_correlation_histogram(python::object& hist,
                              const array<vector<long double>,2>& bins,
                              python::object& ret_bins)
        : _hist(hist), _bins(bins), _ret_bins(ret_bins) {}
115
116
117
118
119
120
121

    template <class Graph, class DegreeSelector1, class DegreeSelector2,
              class WeightMap>
    void operator()(Graph* gp, DegreeSelector1 deg1, DegreeSelector2 deg2,
                    WeightMap weight) const
    {
        Graph& g = *gp;
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
        GetDegreePair put_point;

        typedef typename DegreeSelector1::value_type type1;
        typedef typename DegreeSelector2::value_type type2;

        typedef typename graph_tool::detail::
            select_float_and_larger::apply<type1,type2>::type
            val_type;
        typedef typename property_traits<WeightMap>::value_type count_type;
        typedef Histogram<val_type, count_type, 2> hist_t;

        array<vector<val_type>,2> bins;
        for (size_t i = 0; i < bins.size(); ++i)
        {
            bins[i].resize(_bins[i].size());
            for (size_t j = 0; j < bins[i].size(); ++j)
            {
                // we'll attempt to recover from out of bounds conditions
                try
                {
                    bins[i][j] =
                        numeric_cast<val_type,long double>(_bins[i][j]);
                }
                catch (boost::numeric::negative_overflow&)
                {
                    bins[i][j] = boost::numeric::bounds<val_type>::lowest();
                }
                catch (boost::numeric::positive_overflow&)
                {
                    bins[i][j] = boost::numeric::bounds<val_type>::highest();
                }
            }
            // sort the bins
            sort(bins[i].begin(), bins[i].end());
            // clean bins of zero size
            vector<val_type> temp_bin(1);
            temp_bin[0] = bins[i][0];
            for (size_t j = 1; j < bins[i].size(); ++j)
            {
                if (bins[i][j] > bins[i][j-1])
                    temp_bin.push_back(bins[i][j]);
            }
            bins[i] = temp_bin;
        }

        // find the data range
        pair<type1,type1> range1;
        pair<type2,type2> range2;
        typename graph_traits<Graph>::vertex_iterator vi,vi_end;
        range1.first = boost::numeric::bounds<type1>::highest();
        range1.second = boost::numeric::bounds<type1>::lowest();
        range2.first = boost::numeric::bounds<type2>::highest();
        range2.second = boost::numeric::bounds<type2>::lowest();
        for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
        {
            type1 v1 = deg1(*vi, g);
            type2 v2 = deg2(*vi, g);
            range1.first = min(range1.first, v1);
            range1.second = max(range1.second, v1);
            range2.first = min(range2.first, v2);
            range2.second = max(range2.second, v2);
        }

        boost::array<pair<val_type, val_type>, 2> data_range;
        data_range[0] = range1;
        data_range[1] = range2;

        hist_t hist(bins, data_range);
        SharedHistogram<hist_t> s_hist(hist);
191
192
193
194
195
196
197
198
199

        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;
200
            put_point(v, deg1, deg2, g, weight, s_hist);
201
202
        }
        s_hist.Gather();
203
204
205
206
207
208
209

        bins = hist.GetBins();
        python::list ret_bins;
        ret_bins.append(wrap_vector_owned(bins[0]));
        ret_bins.append(wrap_vector_owned(bins[1]));
        _ret_bins = ret_bins;
        _hist = wrap_multi_array_owned<count_type,2>(hist.GetArray());
210
    }
211
212
213
    python::object& _hist;
    const array<vector<long double>,2>& _bins;
    python::object& _ret_bins;
214
215
216
217
218
219
};


} // graph_tool namespace

#endif