graph_blockmodel.hh 91.3 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-2015 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
22
//
// 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_BLOCKMODEL_HH
#define GRAPH_BLOCKMODEL_HH

#include <cmath>
#include <iostream>
23
#include <queue>
24
#include <boost/math/special_functions/zeta.hpp>
25
#include <boost/functional/hash.hpp>
26

27
#include "config.h"
Tiago Peixoto's avatar
Tiago Peixoto committed
28
29
30
#include <unordered_set>
#include <unordered_map>
#include <tuple>
31

32
#ifdef HAVE_SPARSEHASH
33
34
#include SPARSEHASH_INCLUDE(dense_hash_set)
#include SPARSEHASH_INCLUDE(dense_hash_map)
35
#endif
36

37
#include "../generation/sampler.hh"
38
#include "../generation/dynamic_sampler.hh"
39

40
41
42
43
44
45
46
#ifdef USING_OPENMP
#include <omp.h>
#endif


double spence(double);

47
48
49
namespace graph_tool
{

50
#ifdef HAVE_SPARSEHASH
51
using google::dense_hash_set;
52
using google::dense_hash_map;
53
#endif
54
55
56

using namespace boost;

57
58
59
60
61
62
63
64
65
66
67
68
69
70
template <class Key>
class IdentityArrayPropertyMap
    : public boost::put_get_helper<Key, IdentityArrayPropertyMap<Key>>
{
public:
    typedef std::array<Key, 1> value_type;
    typedef value_type reference;
    typedef Key key_type;
    typedef boost::readable_property_map_tag category;

    inline __attribute__((always_inline))
    const value_type operator[](const key_type& c) const { return {c}; }
};

71
72
73
74
// ====================
// Entropy calculation
// ====================

75
76
77
78
79
80
// Repeated computation of x*log(x) and log(x) actually adds up to a lot of
// time. A significant speedup can be made by caching pre-computed values. This
// is doable since the values of mrse are bounded in [0, 2E], where E is the
// total number of edges in the network. Here we simply grow the cache as
// needed.

81
82
83
84
extern vector<double> __safelog_cache;
extern vector<double> __xlogx_cache;
extern vector<double> __lgamma_cache;

85
86
87
88
89
90
91
92
93
94
template <class Type>
__attribute__((always_inline))
inline double safelog(Type x)
{
    if (x == 0)
        return 0;
    return log(x);
}

__attribute__((always_inline))
95
96
97
98
99
inline double safelog(size_t x)
{
    assert(x < __safelog_cache.size());
    return __safelog_cache[x];
}
100
101

__attribute__((always_inline))
102
103
104
inline double xlogx(size_t x)
{
    if (x >= __xlogx_cache.size())
105
        cout << x << " " << __xlogx_cache.size() << endl;
106
107
108
    assert(x < __xlogx_cache.size());
    return __xlogx_cache[x];
}
109
110

__attribute__((always_inline))
111
112
113
114
115
116
117
inline double lgamma_fast(size_t x)
{
    if (x >= __lgamma_cache.size())
        return lgamma(x);
    return __lgamma_cache[x];
}

118
// polylogarithm and degree-distribution description length (xi)
119
120
121
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

template <class Type>
Type polylog(int n, Type z, Type epsilon=1e-6)
{
    if (n == 2)
        return spence(1 - z);

    int k = 1;
    Type S = 0;
    Type delta = epsilon + 1;
    Type zk = z;
    while (delta > epsilon)
    {
        Type dS = zk / pow(k, n);
        k++;
        zk *= z;
        S += dS;
        delta = dS;
    }
    return S;
}

template <class NType, class EType>
void get_mu_l(NType N, EType E, double& mu, double& l,
              double epsilon=1e-8)
{
    mu = sqrt(polylog<double>(2, 1.) / double(E));
    l = 1. - exp(-double(N) * mu);

    double delta = epsilon + 1;
    while (delta > epsilon)
    {
        double nmu = sqrt(polylog<double>(2, l) / double(E));
        double nl = 1. - exp(-double(N) * mu);

        delta = abs(nmu - mu) + abs(nl - l);
        mu = nmu;
        l = nl;
    }

    l = -log(l);
}

template <class NType, class EType>
double get_xi(NType N, EType E, double epsilon=1e-8)
{
    if (E == 0 || N == 0)
        return 0;
167

168
    double mu = 0, l = 0;
169
    get_mu_l(N, E, mu, l, epsilon);
170
171
172
173
174
175
176
177
178
179
180
181
    double S = double(N) * l + 2 * double(E) * mu;
    return S;
}

template <class NType, class EType>
double get_xi_fast(NType N, EType E)
{
    if (E == 0 || N == 0)
        return 0;
    static const double z2 = boost::math::zeta(2.);
    return 2 * sqrt(z2 * E);
}
182

183
184
// Sparse entropy terms
// ====================
185
186

// "edge" term of the entropy
187
template <class Graph>
188
__attribute__((always_inline))
189
inline double eterm(size_t r, size_t s, size_t mrs, const Graph&)
190
{
191
192
    if (!is_directed::apply<Graph>::type::value && r == s)
        mrs *= 2;
193

194
    double val = xlogx(mrs);
195
196
197
198
199
200
201
202

    if (is_directed::apply<Graph>::type::value || r != s)
        return -val;
    else
        return -val / 2;
}

// "vertex" term of the entropy
203
204
template <class Graph>
inline double vterm(size_t mrp, size_t mrm, size_t wr, bool deg_corr,
205
                    Graph&)
206
207
208
209
210
211
212
{
    double one = 0.5;

    if (is_directed::apply<Graph>::type::value)
        one = 1;

    if (deg_corr)
213
        return one * (xlogx(mrm) + xlogx(mrp));
214
    else
215
        return one * (mrm * safelog(wr) + mrp * safelog(wr));
216
217
218
219
220
}

struct entropy
{
    template <class Graph, class Eprop, class Vprop>
221
222
223
224
225
226
227
    void operator()(Eprop& mrs, Vprop& mrp, Vprop& mrm, Vprop& wr,
                    bool deg_corr, Graph& g, double& S) const
    {
        S = 0;
        for (auto e : edges_range(g))
            S += eterm(source(e, g), target(e, g), mrs[e], g);
        for (auto v : vertices_range(g))
228
            S += vterm(mrp[v], mrm[v], wr[v], deg_corr, g);
229
230
231
    }
};

232
233
// Parallel edges
// ==============
234
235

template <class Vertex, class List, class Graph, class GetNode>
236
double get_parallel_neighbours_entropy(Vertex v, List& us, Graph&,
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
                                       const GetNode& get_node)
{
    double S = 0;
    for (auto& uc : us)
    {
        auto& u = uc.first;
        auto& m = uc.second;
        if (m > 1)
        {
            if (get_node(u) == size_t(v) && !is_directed::apply<Graph>::type::value)
            {
                assert(m % 2 == 0);
                S += lgamma_fast(m/2 + 1);
            }
            else
            {
                S += lgamma_fast(m + 1);
            }
        }
    }
    return S;
}

struct entropy_parallel_edges
{
    template <class Graph, class Weight>
    void operator()(Graph& g, Weight weight, double& S) const
264
265
    {
        S = 0;
266
267
268
269
270
271
272
273
274
275
276
277
278
279
        auto get_node = [](size_t i) {return i;};
        for (auto v : vertices_range(g))
        {
            unordered_map<decltype(v), int> us;
            for (auto e : out_edges_range(v, g))
            {
                auto u = target(e, g);
                if (u < v && !is_directed::apply<Graph>::type::value)
                    continue;
                us[u] += weight[e];
            }

            S += get_parallel_neighbours_entropy(v, us, g, get_node);
        }
280
281
282
283
    }
};


284
// Dense entropy
285
// =============
286

287
288
__attribute__((always_inline))
inline double lbinom(double N, double k)
289
{
290
291
    if (N == 0 || k == 0 || k > N)
        return 0;
292
    return lgamma(N + 1) - lgamma(N - k + 1) - lgamma(k + 1);
293
294
}

295
296
__attribute__((always_inline))
inline double lbinom_fast(int N, int k)
297
{
298
299
    if (N == 0 || k == 0 || k > N)
        return 0;
300
301
    return lgamma_fast(N + 1) - lgamma_fast(N - k + 1) - lgamma_fast(k + 1);
}
302
303


304
305
306
// "edge" term of the entropy
template <class Graph>
__attribute__((always_inline))
307
inline double eterm_dense(size_t r, size_t s, int ers, double wr_r,
308
                          double wr_s, bool multigraph, const Graph&)
309
{
310
311
    // we should not use integers here, since they may overflow
    double nrns;
312

313
314
    if (ers == 0)
        return 0.;
315

316
317
318
    if (r != s || is_directed::apply<Graph>::type::value)
    {
        nrns = wr_r * wr_s;
319
    }
320
    else
321
    {
322
323
324
325
        if (multigraph)
            nrns = (wr_r * (wr_r + 1)) / 2;
        else
            nrns = (wr_r * (wr_r - 1)) / 2;
326
327
    }

328
329
    double S;
    if (multigraph)
330
        S = lbinom(nrns + ers - 1, ers); // do not use lbinom_fast!
331
    else
332
        S = lbinom(nrns, ers);
333
334
    return S;
}
335

336
337
struct entropy_dense
{
338
    template <class Graph, class Eprop, class Vprop>
339
    void operator()(Eprop mrs, Vprop& wr, bool multigraph, Graph& g, double& S) const
340
    {
341
        S = 0;
342
        for (auto e : edges_range(g))
343
        {
344
345
346
            auto r = source(e, g);
            auto s = target(e, g);
            S += eterm_dense(r, s, mrs[e], wr[r], wr[s], multigraph, g);
347
348
349
350
        }
    }
};

351
// ===============
352
// Partition stats
353
// ===============
354
355
356
357
358
359

class partition_stats_t
{
public:

#ifdef HAVE_SPARSEHASH
360
    typedef dense_hash_map<pair<size_t,size_t>, int, std::hash<pair<size_t,size_t>>> map_t;
361
#else
362
    typedef unordered_map<pair<size_t,size_t>, int> map_t;
363
364
365
366
367
#endif

    partition_stats_t() : _enabled(false) {}

    template <class Graph, class Vprop, class Eprop>
368
369
370
371
    partition_stats_t(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
                      bool edges_dl)
        : _enabled(true), _N(N), _E(0), _B(B), _hist(B), _total(B), _ep(B),
          _em(B), _edges_dl(edges_dl)
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
    {

#ifdef HAVE_SPARSEHASH
        for (size_t r = 0; r < B; ++r)
        {
            _hist[r].set_empty_key(make_pair(numeric_limits<size_t>::max(),
                                             numeric_limits<size_t>::max()));
            _hist[r].set_deleted_key(make_pair(numeric_limits<size_t>::max() - 1,
                                               numeric_limits<size_t>::max() - 1));
        }
#endif

        for (auto v : vertices_range(g))
        {
            auto r = b[v];
387
388
389
            size_t kin = in_degreeS()(v, g, eweight);
            size_t kout = out_degreeS()(v, g, eweight);
            _hist[r][make_pair(kin, kout)]++;
390
            _total[r]++;
391
392
393
            _em[r] += kin;
            _ep[r] += kout;
            _E += kout;
394
        }
395
396
397
398
399
400
401
402

        if (!is_directed::apply<Graph>::type::value)
            _E /= 2;

        _actual_B = 0;
        for (size_t r = 0; r < B; ++r)
            if (_total[r] > 0)
                _actual_B++;
403
404
405
406
407
    }

    double get_partition_dl()
    {
        double S = 0;
408
        S += lbinom(_actual_B + _N - 1, _N);
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
        S += lgamma(_N + 1);
        for (auto nr : _total)
            S -= lgamma(nr + 1);
        return S;
    }

    double get_deg_dl(bool ent, bool dl_alt, bool xi_fast)
    {
        double S = 0;
        for (size_t r = 0; r < _B; ++r)
        {
            if (ent)
            {
                for (auto& k_c : _hist[r])
                {
                    double p = k_c.second / double(_total[r]);
                    S -= p * log(p) * _total[r];
                }
            }
            else
            {
                double S1 = 0;

                if (xi_fast)
                {
                    S1 += get_xi_fast(_total[r], _ep[r]);
                    S1 += get_xi_fast(_total[r], _em[r]);
                }
                else
                {
                    S1 += get_xi(_total[r], _ep[r]);
                    S1 += get_xi(_total[r], _em[r]);
                }

                S1 += lgamma(_total[r] + 1);
                for (auto& k_c : _hist[r])
                    S1 -= lgamma(k_c.second + 1);

                if (dl_alt)
                {
                    double S2 = 0;
                    S2 += lbinom(_total[r] + _ep[r] - 1, _ep[r]);
                    S2 += lbinom(_total[r] + _em[r] - 1, _em[r]);
                    S += min(S1, S2);
                }
                else
                {
                    S += S1;
                }
            }
        }
        return S;
    }

    template <class Graph, class OStats>
464
    double get_delta_dl(size_t, size_t r, size_t nr, OStats&, Graph&)
465
466
467
468
469
470
471
472
    {
        if (r == nr)
            return 0;

        double S_b = 0, S_a = 0;
        S_b += -lgamma_fast(_total[r] + 1) - lgamma_fast(_total[nr] + 1);
        S_a += -lgamma_fast(_total[r]    ) - lgamma_fast(_total[nr] + 2);

473
474
475
476
477
478
479
        int dB = 0;
        if (_total[r] == 1)
            dB--;
        if (_total[nr] == 0)
            dB++;

        if (dB != 0)
480
        {
481
482
            S_b += lbinom(_actual_B + _N - 1, _N);
            S_a += lbinom(_actual_B + dB + _N - 1, _N);
483

484
485
486
            if (_edges_dl)
            {
                auto get_x = [](size_t B) -> size_t
487
                {
488
489
490
491
                    if (is_directed::apply<Graph>::type::value)
                        return B * B;
                    else
                        return (B * (B + 1)) / 2;
492
493
                };

494
495
496
497
                S_b += lbinom(get_x(_actual_B) + _E - 1, _E);
                S_a += lbinom(get_x(_actual_B + dB) + _E - 1, _E);
            }
        }
498

499
500
501
502
503
504
505
506
507
        return S_a - S_b;
    }

    template <class Graph, class EWeight, class OStats>
    double get_delta_deg_dl(size_t v, size_t r, size_t nr, EWeight& eweight,
                            OStats&, Graph& g)
    {
        if (r == nr)
            return 0;
508

509
        double S_b = 0, S_a = 0;
510

511
512
        int kin = in_degreeS()(v, g, eweight);
        int kout = out_degreeS()(v, g, eweight);
513

514
515
516
517
518
519
520
        auto get_Se = [&](size_t s, int delta, int kin, int kout) -> double
            {
                double S = 0;
                S += get_xi_fast(_total[s] + delta, _em[s] + kin);
                S += get_xi_fast(_total[s] + delta, _ep[s] + kout);
                return S;
            };
521

522
523
        S_b += get_Se(r,  0,    0,     0) + get_Se(nr, 0,   0,    0);
        S_a += get_Se(r, -1, -kin, -kout) + get_Se(nr, 1, kin, kout);
524

525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
        auto get_Sr = [&](size_t s, int delta) -> double
            {
                return lgamma_fast(_total[s] + delta + 1);
            };

        S_b += get_Sr(r,  0) + get_Sr(nr, 0);
        S_a += get_Sr(r, -1) + get_Sr(nr, 1);


        auto get_Sk = [&](size_t s, pair<size_t, size_t>& deg, int delta) -> double
            {
                size_t nd = 0;
                auto iter = _hist[s].find(deg);
                if (iter != _hist[s].end())
                    nd = iter->second;

                return -lgamma_fast(nd + delta + 1);
            };

        auto deg = make_pair(size_t(kin), size_t(kout));
        S_b += get_Sk(r, deg,  0) + get_Sk(nr, deg, 0);
        S_a += get_Sk(r, deg, -1) + get_Sk(nr, deg, 1);
547
548
549
550
551

        return S_a - S_b;
    }

    template <class Graph, class OStats>
552
553
    void move_vertex(size_t v, size_t r, size_t nr, bool deg_corr, OStats&,
                     Graph& g, size_t kin = 0, size_t kout = 0)
554
555
556
557
558
559
560
    {
        if (r == nr)
            return;

        _total[r]--;
        _total[nr]++;

561
562
563
564
565
        if (_total[r] == 0)
            _actual_B--;
        if (_total[nr] == 1)
            _actual_B++;

566
567
        if (deg_corr)
        {
568
569
570
571
572
573
            if (kin + kout == 0)
            {
                kin = in_degreeS()(v, g);
                kout = out_degreeS()(v, g);
            }
            auto deg = make_pair(kin, kout);
574
575
576
577
578
579
580
581
582
583
584
585
586
            auto iter = _hist[r].find(deg);
            iter->second--;
            if (iter->second == 0)
                _hist[r].erase(iter);
            _hist[nr][deg]++;
            _em[r] -= deg.first;
            _ep[r] -= deg.second;
            _em[nr] += deg.first;
            _ep[nr] += deg.second;
        }
    }

    bool is_enabled() { return _enabled; }
587
    void set_enabled(bool enabled) { _enabled = enabled; }
588
589
590
591

private:
    bool _enabled;
    size_t _N;
592
    size_t _E;
593
    size_t _B;
594
    size_t _actual_B;
595
596
597
598
    vector<map_t> _hist;
    vector<int> _total;
    vector<int> _ep;
    vector<int> _em;
599
    bool _edges_dl;
600
};
601

602
603
604
605
// ===============================
// Block moves
// ===============================

606
607
// this structure speeds up the access to the edges between given blocks,
// since we're using an adjacency list to store the block structure (the emat_t
608
609
610
611
612
613
614
615
616
617
618
619
620
621
// is simply a corresponding adjacency matrix)
struct get_emat_t
{
    template <class Graph>
    struct apply
    {
        typedef multi_array<pair<typename graph_traits<Graph>::edge_descriptor, bool>, 2> type;
    };
};


struct create_emat
{
    template <class Graph>
622
    void operator()(Graph& g, boost::any& oemap) const
623
624
    {
        typedef typename get_emat_t::apply<Graph>::type emat_t;
625
        size_t B = num_vertices(g);
626
627
628
629
630
        emat_t emat(boost::extents[B][B]);

        for (size_t i = 0; i < B; ++i)
            for (size_t j = 0; j < B; ++j)
                emat[i][j].second = false;
631
        for (auto e : edges_range(g))
632
        {
633
            if (source(e, g) >= B || target(e, g) >= B)
634
                throw GraphException("incorrect number of blocks when creating emat!");
635
            emat[source(e, g)][target(e, g)] = make_pair(e, true);
636
            if (!is_directed::apply<Graph>::type::value)
637
                emat[target(e, g)][source(e, g)] = make_pair(e, true);
638
639
640
641
642
643
        }

        oemap = emat;
    }
};

644
645
646
647
648
template <class Graph>
inline __attribute__((always_inline))
pair<typename graph_traits<Graph>::edge_descriptor, bool>
get_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
649
       const typename get_emat_t::apply<Graph>::type& emat, const Graph&)
650
651
652
653
654
655
656
657
658
659
{
    return emat[r][s];
}

template <class Graph>
inline __attribute__((always_inline))
void
put_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename graph_traits<Graph>::edge_descriptor& e,
660
       typename get_emat_t::apply<Graph>::type& emat, const Graph&)
661
662
663
664
665
666
667
668
669
670
671
{
    emat[r][s] = make_pair(e, true);
    if (!is_directed::apply<Graph>::type::value && r != s)
        emat[s][r] = make_pair(e, true);
}

template <class Graph>
inline __attribute__((always_inline))
void
remove_me(typename graph_traits<Graph>::vertex_descriptor r,
          typename graph_traits<Graph>::vertex_descriptor s,
672
673
          const typename graph_traits<Graph>::edge_descriptor&,
          typename get_emat_t::apply<Graph>::type& emat, Graph&,
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
          bool delete_edge=true)
{
    if (!delete_edge)
    {
        emat[r][s].second = false;
        if (!is_directed::apply<Graph>::type::value && r != s)
            emat[s][r].second = false;
    }
}


// this structure speeds up the access to the edges between given blocks, since
// we're using an adjacency list to store the block structure (this is like
// emat_t above, but takes less space and is slower)
struct get_ehash_t
{
    template <class Graph>
    struct apply
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
#ifdef HAVE_SPARSEHASH
696
        typedef dense_hash_map<vertex_t, edge_t, std::hash<vertex_t>> map_t;
697
698
699
700
701
702
703
704
705
706
707
708
709
#else
        typedef unordered_map<vertex_t, edge_t> map_t;
#endif
        typedef vector<map_t> type;
    };
};


template<class Graph>
inline __attribute__((always_inline))
pair<typename graph_traits<Graph>::edge_descriptor, bool>
get_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
710
       const typename get_ehash_t::apply<Graph>::type& ehash, const Graph&)
711
{
712
    assert(r < ehash.size());
713
714
    const auto& map = ehash[r];
    auto iter = map.find(s);
715
716
717
718
719
720
721
722
723
724
725
726
    if (iter == map.end())
        return (make_pair(typename graph_traits<Graph>::edge_descriptor(), false));
    return make_pair(iter->second, true);
}

template<class Graph>
inline __attribute__((always_inline))
void
put_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename graph_traits<Graph>::edge_descriptor& e,
       typename get_ehash_t::apply<Graph>::type& ehash,
727
       const Graph&)
728
{
729
    assert(r < ehash.size());
730
731
732
733
734
735
736
737
738
739
740
741
742
743
    ehash[r][s] = e;
    if (!is_directed::apply<Graph>::type::value)
        ehash[s][r] = e;
}

template<class Graph>
inline __attribute__((always_inline))
void
remove_me(typename graph_traits<Graph>::vertex_descriptor r,
          typename graph_traits<Graph>::vertex_descriptor s,
          const typename graph_traits<Graph>::edge_descriptor& e,
          typename get_ehash_t::apply<Graph>::type& ehash, Graph& bg,
          bool delete_edge=true)
{
744
    assert(r < ehash.size());
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    ehash[r].erase(s);
    if (!is_directed::apply<Graph>::type::value)
        ehash[s].erase(r);
    if (delete_edge)
        remove_edge(e, bg);
}

struct create_ehash
{
    template <class Graph>
    void operator()(Graph& g, boost::any& oemap) const
    {
        typedef typename get_ehash_t::apply<Graph>::type emat_t;
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

        emat_t emat(num_vertices(g));

#ifdef HAVE_SPARSEHASH
763
        for (auto v : vertices_range(g))
764
        {
765
            emat[v].set_empty_key(numeric_limits<vertex_t>::max());
766
            emat[v].set_deleted_key(numeric_limits<vertex_t>::max() - 1);
767
768
769
        }
#endif

770
771
        for (auto e : edges_range(g))
            put_me(source(e, g), target(e, g), e, emat, g);
772

773
774
775
776
#ifdef HAVE_SPARSEHASH
        for (auto v : vertices_range(g))
            emat[v].resize(0);
#endif
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
        oemap = emat;
    }
};

template <class Vertex, class Eprop, class Emat, class BGraph>
__attribute__((always_inline))
inline size_t get_mrs(Vertex r, Vertex s, const Eprop& mrs, Emat& emat,
                      BGraph& bg)
{
    const pair<typename graph_traits<BGraph>::edge_descriptor, bool> me =
        get_me(r, s, emat, bg);
    if (me.second)
        return mrs[me.first];
    else
        return 0;
}

794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
struct standard_neighbours_policy
{
    template <class Graph, class Vertex>
    IterRange<typename out_edge_iteratorS<Graph>::type>
    get_out_edges(Vertex v, Graph& g) const
    {
        return out_edges_range(v, g);
    }

    template <class Graph, class Vertex>
    IterRange<typename in_edge_iteratorS<Graph>::type>
    get_in_edges(Vertex v, Graph& g) const
    {
        return in_edges_range(v, g);
    }

    template <class Graph, class Vertex, class Weight>
    int get_out_degree(Vertex& v, Graph& g, Weight& eweight) const
    {
        return out_degreeS()(v, g, eweight);
    }

    template <class Graph, class Vertex, class Weight>
    int get_in_degree(Vertex& v, Graph& g, Weight& eweight) const
    {
        return in_degreeS()(v, g, eweight);
    }
};

823
824
// remove a vertex from its current block
template <class Graph, class BGraph, class Eprop, class Vprop, class EWprop,
825
826
          class VWprop, class EMat, class OStats,
          class NPolicy = standard_neighbours_policy>
827
828
void remove_vertex(size_t v, Eprop& mrs, Vprop& mrp, Vprop& mrm, Vprop& wr,
                   Vprop& b, const EWprop& eweight, const VWprop& vweight,
829
                   Graph& g, BGraph& bg, EMat& emat, OStats& overlap_stats,
830
                   const NPolicy& npolicy = NPolicy())
831
832
833
834
{
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    vertex_t r = b[v];

835
    int self_weight = 0;
836
    for (auto e : npolicy.get_out_edges(v, g))
837
    {
838
        vertex_t u = target(e, g);
839
840
        vertex_t s = b[u];

841
842
843
        // if (!emat[r][s].second)
        //     throw GraphException("no edge? " + lexical_cast<string>(r) +
        //                          " " + lexical_cast<string>(s));
844

845
        auto me = get_me(r, s, emat, bg).first;
846

847
        size_t ew = eweight[e];
848
849
850
851
852
853
854
855
856
857
858
859
        if (u == v && !is_directed::apply<Graph>::type::value)
        {
            self_weight += ew;
        }
        else
        {
            mrs[me] -= ew;

            assert(mrs[me] >= 0);

            mrp[r] -= ew;
            mrm[s] -= ew;
860

861
862
863
864
            if (mrs[me] == 0)
                remove_me(r, s, me, emat, bg);
        }
    }
865

866
867
868
869
870
871
872
873
    if (self_weight > 0)
    {
        assert(self_weight % 2 == 0);
        auto me = get_me(r, r, emat, bg).first;
        mrs[me] -= self_weight / 2;
        mrp[r] -= self_weight / 2;
        mrm[r] -= self_weight / 2;
        assert(mrs[me] >= 0);
874
        if (mrs[me] == 0)
875
            remove_me(r, r, me, emat, bg);
876
877
    }

878
    for (auto e : npolicy.get_in_edges(v, g))
879
    {
880
881
882
883
        vertex_t u = source(e, g);
        if (u == v)
            continue;
        vertex_t s = b[u];
884

885
886
887
        // if (!emat[s][r].second)
        //     throw GraphException("no edge? " + lexical_cast<string>(s) +
        //                          " " + lexical_cast<string>(r));
888

889
890
        typename graph_traits<BGraph>::edge_descriptor me =
            get_me(s, r, emat, bg).first;
891

892
893
        size_t ew = eweight[e];
        mrs[me] -= ew;
894

895
896
        mrp[s] -= ew;
        mrm[r] -= ew;
897

898
899
        if (mrs[me] == 0)
            remove_me(s, r, me, emat, bg);
900
901
    }

902
903
904
905
906
907
    if (!overlap_stats.is_enabled())
    {
        wr[r] -= vweight[v];
    }
    else
    {
908
        overlap_stats.remove_half_edge(v, r, b, g);
909
910
        wr[r] = overlap_stats.get_block_size(r);
    }
911
912
913
914
}

// add a vertex to block rr
template <class Graph, class BGraph, class Eprop, class Vprop, class EWprop,
915
916
          class VWprop, class EMat, class OStats,
          class NPolicy = standard_neighbours_policy>
917
void add_vertex(size_t v, size_t r, Eprop& mrs, Vprop& mrp, Vprop& mrm,
918
                Vprop& wr, Vprop& b, const EWprop& eweight,
919
920
                const VWprop& vweight, Graph& g, BGraph& bg, EMat& emat,
                OStats& overlap_stats, const NPolicy& npolicy = NPolicy())
921
922
923
{
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

924
    int self_weight = 0;
925
    for (auto e : npolicy.get_out_edges(v, g))
926
    {
927
        vertex_t u = target(e, g);
928
929
930
931
932
933
934
935
936
        vertex_t s;

        if (u != v)
            s = b[u];
        else
            s = r;

        typename graph_traits<BGraph>::edge_descriptor me;

937
        auto mep = get_me(r, s, emat, bg);
938
939
940

        if (!mep.second)
        {
941
942
943
            mep = add_edge(r, s, bg);
            put_me(r, s, mep.first, emat, bg);
            mrs[mep.first] = 0;
944
        }
945
        me = mep.first;
946

947
        size_t ew = eweight[e];
948

949
950
951
952
953
954
955
956
957
958
959
        if (u == v && !is_directed::apply<Graph>::type::value)
        {
            self_weight += ew;
        }
        else
        {
            mrs[me] += ew;
            mrp[r] += ew;
            mrm[s] += ew;
        }
    }
960

961
962
963
964
965
966
967
968
    if (self_weight > 0)
    {
        assert(self_weight % 2 == 0);
        auto me = get_me(r, r, emat, bg).first;
        mrs[me] += self_weight / 2;
        mrp[r] += self_weight / 2;
        mrm[r] += self_weight / 2;
        assert(mrs[me] >= 0);
969
970
    }

971
    for (auto e : npolicy.get_in_edges(v, g))
972
    {
973
        vertex_t u = source(e, g);
974
975
976
977
978
979
        if (u == v)
            continue;

        vertex_t s = b[u];

        typename graph_traits<BGraph>::edge_descriptor me;
980
981
982
        pair<typename graph_traits<BGraph>::edge_descriptor, bool> mep =
                get_me(s, r, emat, bg);

983
984
985

        if (!mep.second)
        {
986
987
988
            mep = add_edge(s, r, bg);
            put_me(s, r, mep.first, emat, bg);
            mrs[mep.first] = 0;
989
        }
990
        me = mep.first;
991

992
        size_t ew = eweight[e];
993
994
995
996
997
998
999

        mrs[me] += ew;

        mrp[s] += ew;
        mrm[r] += ew;
    }

1000
    if (!overlap_stats.is_enabled())
For faster browsing, not all history is shown. View entire blame