graph_blockmodel_util.hh 42.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.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_BLOCKMODEL_UTIL_HH
#define GRAPH_BLOCKMODEL_UTIL_HH

#include "config.h"

#include <tuple>

#include "hash_map_wrap.hh"

#include "../generation/sampler.hh"
#include "../generation/dynamic_sampler.hh"
29
#include "graph_neighbour_sampler.hh"
30
31

#include "util.hh"
32
#include "int_part.hh"
33
34
35
#include "cache.hh"

#include <boost/multi_array.hpp>
36
#include <boost/math/special_functions/gamma.hpp>
37
38
39
40
41
42
43
44

#ifdef USING_OPENMP
#include <omp.h>
#endif

namespace graph_tool
{

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
// tuple utils
template <size_t i, class T>
void add_to_tuple_imp(T&)
{
}

template <size_t i, class T, class Ti, class... Ts>
void add_to_tuple_imp(T& tuple, Ti&& v, Ts&&... vals)
{
    get<i>(tuple) += v;
    add_to_tuple_imp<i+1>(tuple, vals...);
}

template <class T, class... Ts>
void add_to_tuple(T& tuple, Ts&&... vals)
{
    add_to_tuple_imp<0>(tuple, vals...);
}

namespace detail {
   template <class F, class Tuple, std::size_t... I>
   constexpr decltype(auto) tuple_apply_impl(F&& f, Tuple&& t,
                                             std::index_sequence<I...>)
   {
       return f(std::get<I>(std::forward<Tuple>(t))...);
   }
} // namespace detail

template <class F, class Tuple>
constexpr decltype(auto) tuple_apply(F&& f, Tuple&& t)
{
    return detail::tuple_apply_impl
        (std::forward<F>(f), std::forward<Tuple>(t),
         std::make_index_sequence<std::tuple_size<std::decay_t<Tuple>>{}>{});
}

81
82
83
84
// ====================
// Entropy calculation
// ====================

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
enum deg_dl_kind
{
    ENT,
    UNIFORM,
    DIST
};

struct entropy_args_t
{
    bool dense;
    bool multigraph;
    bool exact;
    bool adjacency;
    bool partition_dl;
    bool degree_dl;
    deg_dl_kind  degree_dl_kind;
    bool edges_dl;
};
103
104
105
106

// Sparse entropy terms
// ====================

107
108
109
110
111
112
113
114
115
116
117
118
// exact microcanonical deg-corr entropy
template <class Graph>
inline double eterm_exact(size_t r, size_t s, size_t mrs, const Graph&)
{
    double val = lgamma_fast(mrs + 1);

    if (is_directed::apply<Graph>::type::value || r != s)
    {
        return -val;
    }
    else
    {
119
#ifndef __clang__
120
        constexpr double log_2 = log(2);
121
122
123
#else
        const double log_2 = log(2);
#endif
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
        return -val - mrs * log_2;
    }
}

template <class Graph>
inline double vterm_exact(size_t mrp, size_t mrm, size_t wr, bool deg_corr,
                          const Graph&)
{
    if (deg_corr)
    {
        if (is_directed::apply<Graph>::type::value)
            return lgamma_fast(mrp + 1) + lgamma_fast(mrm + 1);
        else
            return lgamma_fast(mrp + 1);
    }
    else
    {
        if (is_directed::apply<Graph>::type::value)
            return (mrp + mrm) * safelog(wr);
        else
            return mrp * safelog(wr);
    }
}

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
// "edge" term of the entropy
template <class Graph>
inline double eterm(size_t r, size_t s, size_t mrs, const Graph&)
{
    if (!is_directed::apply<Graph>::type::value && r == s)
        mrs *= 2;

    double val = xlogx(mrs);

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

// "vertex" term of the entropy
template <class Graph>
inline double vterm(size_t mrp, size_t mrm, size_t wr, bool deg_corr,
                    Graph&)
{
    double one = 0.5;

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

    if (deg_corr)
        return one * (xlogx(mrm) + xlogx(mrp));
    else
        return one * (mrm * safelog(wr) + mrp * safelog(wr));
}

179
180


181
182
183
184
185
186
187
188
189
190
191
192
193
194
// Dense entropy
// =============

// "edge" term of the entropy
template <class Graph>
inline double eterm_dense(size_t r, size_t s, int ers, double wr_r,
                          double wr_s, bool multigraph, const Graph&)
{
    // we should not use integers here, since they may overflow
    double nrns;

    if (ers == 0)
        return 0.;

195
196
    assert(wr_r + wr_s > 0);

197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
    if (r != s || is_directed::apply<Graph>::type::value)
    {
        nrns = wr_r * wr_s;
    }
    else
    {
        if (multigraph)
            nrns = (wr_r * (wr_r + 1)) / 2;
        else
            nrns = (wr_r * (wr_r - 1)) / 2;
    }

    double S;
    if (multigraph)
        S = lbinom(nrns + ers - 1, ers); // do not use lbinom_fast!
    else
        S = lbinom(nrns, ers);
    return S;
}

217
218
// Weighted entropy terms

219
// exponential
220
221
222
223
224
225
226
227
228
template <class DT>
double positive_w_log_P(DT N, double x, double alpha, double beta)
{
    if (N == 0)
        return 0.;
    return boost::math::lgamma(N + alpha) - boost::math::lgamma(alpha)
        + alpha * log(beta) - (alpha + N) * log(beta + x);
}

229
// normal
230
231
232
233
234
235
236
237
238
239
240
241
242
243
template <class DT>
double signed_w_log_P(DT N, double x, double v, double m0, double k0, double v0,
                      double nu0)
{
    if (N == 0)
        return 0.;
    auto k_n = k0 + N;
    auto nu_n = nu0 + N;
    auto v_n = (v0 * nu0 + v + ((N * k0)/(k0 + N)) * pow(m0 - x/N, 2.)) / nu_n;
    return boost::math::lgamma(nu_n/2.) - boost::math::lgamma(nu0/2.)
        + (log(k0) - log(k_n))/2. + (nu0 / 2.) * log(nu0 * v0)
        - (nu_n / 2.) * log(nu_n * v_n) - (N/2.) * log(M_PI);
}

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
// discrete: geometric
template <class DT>
double geometric_w_log_P(DT N, double x, double alpha, double beta)
{
    if (N == 0)
        return 0.;
    return boost::math::lgamma(alpha + beta)
        + boost::math::lgamma(alpha + 1) + boost::math::lgamma(beta + 1)
        - boost::math::lgamma(alpha) - boost::math::lgamma(beta)
        - boost::math::lgamma(x + alpha + beta + 1);
}

// discrete: Poisson
template <class DT>
double poisson_w_log_P(DT N, double x, double alpha, double beta)
{
    if (N == 0)
        return 0.;
    return boost::math::lgamma(x + alpha)
        - boost::math::lgamma(x + 1) - boost::math::lgamma(alpha)
        + alpha * log(beta) - (x + alpha) * log(beta + 1);
}
266

267
268
269
270
// ===============
// Partition stats
// ===============

271
272
constexpr size_t null_group = std::numeric_limits<size_t>::max();

273
274
275
276
277
278
279
280
281
282
typedef vprop_map_t<std::vector<std::tuple<size_t, size_t, size_t>>>::type
    degs_map_t;

struct simple_degs_t {};

template <class Graph, class Vprop, class Eprop>
auto get_degs(size_t v, Vprop& vweight, Eprop& eweight, const simple_degs_t&,
              Graph& g)
{
    std::array<std::tuple<size_t, size_t, size_t>, 1>
283
284
285
        degs {{make_tuple(size_t(in_degreeS()(v, g, eweight)),
                          size_t(out_degreeS()(v, g, eweight)),
                          size_t(vweight[v]))}};
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
    return degs;
}

template <class Graph, class Vprop, class Eprop>
auto& get_degs(size_t v, Vprop&, Eprop&, typename degs_map_t::unchecked_t& degs,
               Graph&)
{
    return degs[v];
}

class partition_stats_t
{
public:

    typedef gt_hash_map<pair<size_t,size_t>, int> map_t;

    template <class Graph, class Vprop, class VWprop, class Eprop, class Degs,
              class Mprop, class Vlist>
    partition_stats_t(Graph& g, Vprop& b, Vlist& vlist, size_t E, size_t B,
                      VWprop& vweight, Eprop& eweight, Degs& degs,
306
307
                      const Mprop& ignore_degree, std::vector<size_t>& bmap,
                      bool allow_empty)
308
        : _bmap(bmap), _N(0), _E(E), _total_B(B), _allow_empty(allow_empty)
309
310
311
    {
        for (auto v : vlist)
        {
312
313
314
            if (vweight[v] == 0)
                continue;

315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
            auto r = get_r(b[v]);

            auto&& ks = get_degs(v, vweight, eweight, degs, g);
            if (v >= _ignore_degree.size())
                _ignore_degree.resize(v + 1, 0);
            _ignore_degree[v] = ignore_degree[v];
            for (auto& k : ks)
            {
                size_t kin = get<0>(k);
                size_t kout = get<1>(k);
                size_t n = get<2>(k);
                if (_ignore_degree[v] == 2)
                    kout = 0;
                if (_ignore_degree[v] != 1)
                {
                    _hist[r][make_pair(kin, kout)] += n;
                    _em[r] += kin * n;
                    _ep[r] += kout * n;
                }
                _total[r] += n;
                _N += n;
            }
        }

        _actual_B = 0;
        for (auto n : _total)
        {
            if (n > 0)
                _actual_B++;
        }
    }

    size_t get_r(size_t r)
    {
        constexpr size_t null =
            std::numeric_limits<size_t>::max();
        if (r >= _bmap.size())
            _bmap.resize(r + 1, null);
        size_t nr = _bmap[r];
        if (nr == null)
            nr = _bmap[r] = _hist.size();
        if (nr >= _hist.size())
        {
            _hist.resize(nr + 1);
            _total.resize(nr + 1);
            _ep.resize(nr + 1);
            _em.resize(nr + 1);
        }
        return nr;
    }

    double get_partition_dl()
    {
        double S = 0;
369
        if (_allow_empty)
370
            S += lbinom(_total_B + _N - 1, _N);
371
372
        else
            S += lbinom(_N - 1, _actual_B - 1);
373
        S += lgamma_fast(_N + 1);
374
        for (auto nr : _total)
375
            S -= lgamma_fast(nr + 1);
376
        S += safelog(_N);
377
378
379
        return S;
    }

380
    double get_deg_dl_ent()
381
382
383
384
    {
        double S = 0;
        for (size_t r = 0; r < _ep.size(); ++r)
        {
385
            size_t total = 0;
386
            for (auto& k_c : _hist[r])
387
            {
388
                S -= xlogx(k_c.second);
389
390
391
                total += k_c.second;
            }
            S += xlogx(total);
392
393
394
        }
        return S;
    }
395

396
397
398
399
400
401
402
403
404
405
    double get_deg_dl_uniform()
    {
        double S = 0;
        for (size_t r = 0; r < _ep.size(); ++r)
        {
            S += lbinom(_total[r] + _ep[r] - 1, _ep[r]);
            S += lbinom(_total[r] + _em[r] - 1, _em[r]);
        }
        return S;
    }
406

407
408
409
410
411
412
413
    double get_deg_dl_dist()
    {
        double S = 0;
        for (size_t r = 0; r < _ep.size(); ++r)
        {
            S += log_q(_ep[r], _total[r]);
            S += log_q(_em[r], _total[r]);
414

415
            size_t total = 0;
416
            for (auto& k_c : _hist[r])
417
            {
418
                S -= lgamma_fast(k_c.second + 1);
419
420
421
                total += k_c.second;
            }
            S += lgamma_fast(total + 1);
422
423
424
425
        }
        return S;
    }

426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
    double get_deg_dl(int kind)
    {
        switch (kind)
        {
        case deg_dl_kind::ENT:
            return get_deg_dl_ent();
        case deg_dl_kind::UNIFORM:
            return get_deg_dl_uniform();
        case deg_dl_kind::DIST:
            return get_deg_dl_dist();
        default:
            return numeric_limits<double>::quiet_NaN();
        }
    }

441
    template <class VProp>
442
    double get_delta_partition_dl(size_t v, size_t r, size_t nr, VProp& vweight)
443
444
445
446
    {
        if (r == nr)
            return 0;

447
448
449
450
451
452
453
        if (r != null_group)
            r = get_r(r);

        if (nr != null_group)
            nr = get_r(nr);

        int n = vweight[v];
454
        if (n == 0)
455
456
457
458
459
460
        {
            if (r == null_group)
                n = 1;
            else
                return 0;
        }
461
462
463

        double S_b = 0, S_a = 0;

464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
        if (r != null_group)
        {
            S_b += -lgamma_fast(_total[r] + 1);
            S_a += -lgamma_fast(_total[r] - n + 1);
        }

        if (nr != null_group)
        {
            S_b += -lgamma_fast(_total[nr] + 1);
            S_a += -lgamma_fast(_total[nr] + n + 1);
        }

        int dN = 0;
        if (r == null_group)
            dN += n;
        if (nr == null_group)
            dN -= n;

        S_b += lgamma_fast(_N + 1);
        S_a += lgamma_fast(_N + dN + 1);
484
485

        int dB = 0;
486
        if (r != null_group && _total[r] == n)
487
            dB--;
488
        if (nr != null_group && _total[nr] == 0)
489
490
            dB++;

491
        if ((dN != 0 || dB != 0) && !_allow_empty)
492
        {
493
494
            S_b += lbinom_fast(_N - 1, _actual_B - 1);
            S_a += lbinom_fast(_N - 1 + dN, _actual_B + dB - 1);
495
496
        }

497
498
499
500
501
502
        if (dN != 0)
        {
            S_b += safelog(_N);
            S_a += safelog(_N + dN);
        }

503
504
505
        return S_a - S_b;
    }

506
    template <class VProp, class Graph>
507
    double get_delta_edges_dl(size_t v, size_t r, size_t nr, VProp& vweight,
508
                              size_t actual_B, Graph&)
509
    {
510
        if (r == nr || _allow_empty)
511
512
            return 0;

513
514
515
516
        if (r != null_group)
            r = get_r(r);
        if (nr != null_group)
            nr = get_r(nr);
517
518
519
520
521
522

        double S_b = 0, S_a = 0;

        int n = vweight[v];

        if (n == 0)
523
524
525
526
527
528
        {
            if (r == null_group)
                n = 1;
            else
                return 0;
        }
529
530

        int dB = 0;
531
        if (r != null_group && _total[r] == n)
532
            dB--;
533
        if (nr != null_group && _total[nr] == 0)
534
535
536
537
538
539
540
541
542
543
544
545
            dB++;

        if (dB != 0)
        {
            auto get_x = [](size_t B)
                {
                    if (is_directed::apply<Graph>::type::value)
                        return B * B;
                    else
                        return (B * (B + 1)) / 2;
                };

546
547
            S_b += lbinom(get_x(actual_B) + _E - 1, _E);
            S_a += lbinom(get_x(actual_B + dB) + _E - 1, _E);
548
549
550
551
552
553
554
        }

        return S_a - S_b;
    }

    template <class Graph, class VProp, class EProp, class Degs>
    double get_delta_deg_dl(size_t v, size_t r, size_t nr, VProp& vweight,
555
                            EProp& eweight, Degs& degs, Graph& g, int kind)
556
    {
557
        if (r == nr || _ignore_degree[v] == 1 || vweight[v] == 0)
558
            return 0;
559
560
561
562
        if (r != null_group)
            r = get_r(r);
        if (nr != null_group)
            nr = get_r(nr);
563
        auto&& ks = get_degs(v, vweight, eweight, degs, g);
564
565
566
567
568
569
570
571
572
        auto* _ks = &ks;
        typename std::remove_reference<decltype(ks)>::type nks;
        if (_ignore_degree[v] == 2)
        {
            nks = ks;
            for (auto& k : nks)
                get<1>(k) = 0;
            _ks = &nks;
        }
573
        double dS = 0;
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
        switch (kind)
        {
        case deg_dl_kind::ENT:
            if (r != null_group)
                dS += get_delta_deg_dl_ent_change(r,  *_ks, -1);
            if (nr != null_group)
                dS += get_delta_deg_dl_ent_change(nr, *_ks, +1);
            break;
        case deg_dl_kind::UNIFORM:
            if (r != null_group)
                dS += get_delta_deg_dl_uniform_change(v, r,  *_ks, -1);
            if (nr != null_group)
                dS += get_delta_deg_dl_uniform_change(v, nr, *_ks, +1);
            break;
        case deg_dl_kind::DIST:
            if (r != null_group)
                dS += get_delta_deg_dl_dist_change(v, r,  *_ks, -1);
            if (nr != null_group)
                dS += get_delta_deg_dl_dist_change(v, nr, *_ks, +1);
            break;
        default:
            dS = numeric_limits<double>::quiet_NaN();
        }
597
598
599
        return dS;
    }

600
601
    template <class Ks>
    double get_delta_deg_dl_ent_change(size_t r, Ks& ks, int diff)
602
    {
603
604
605
606
607
608
609
610
611
612
        int nr = _total[r];
        auto get_Sk = [&](size_t s, pair<size_t, size_t>& deg, int delta)
            {
                int nd = 0;
                auto iter = _hist[s].find(deg);
                if (iter != _hist[s].end())
                    nd = iter->second;
                assert(nd + delta >= 0);
                return -xlogx(nd + delta);
            };
613

614
615
616
617
618
619
620
621
622
623
624
625
626
627
        double S_b = 0, S_a = 0;
        int dn = 0;
        for (auto& k : ks)
        {
            size_t kin = get<0>(k);
            size_t kout = get<1>(k);
            int nk = get<2>(k);
            dn += diff * nk;
            auto deg = make_pair(kin, kout);
            S_b += get_Sk(r, deg,         0);
            S_a += get_Sk(r, deg, diff * nk);
        }
        S_b += xlogx(nr);
        S_a += xlogx(nr + dn);
628

629
630
631
632
633
        return S_a - S_b;
    }

    template <class Ks>
    double get_delta_deg_dl_uniform_change(size_t v, size_t r, Ks& ks, int diff)
634
    {
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
        auto get_Se = [&](int dn, int dkin, int dkout)
            {
                double S = 0;
                S += lbinom(_total[r] + dn + _ep[r] - 1 + dkout, _ep[r] + dkout);
                S += lbinom(_total[r] + dn + _em[r] - 1 + dkin,  _em[r] + dkin);
                return S;
            };

        double S_b = 0, S_a = 0;
        int tkin = 0, tkout = 0, n = 0;
        for (auto& k : ks)
        {
            size_t kin = get<0>(k);
            size_t kout = get<1>(k);
            int nk = get<2>(k);
            tkin += kin * nk;
            if (_ignore_degree[v] != 2)
                tkout += kout * nk;
            n += nk;
        }

        S_b += get_Se(       0,           0,            0);
        S_a += get_Se(diff * n, diff * tkin, diff * tkout);
        return S_a - S_b;
    }
660
661

    template <class Ks>
662
    double get_delta_deg_dl_dist_change(size_t v, size_t r, Ks& ks, int diff)
663
    {
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691

        auto get_Se = [&](int delta, int kin, int kout)
            {
                double S = 0;
                assert(_total[r] + delta >= 0);
                assert(_em[r] + kin >= 0);
                assert(_ep[r] + kout >= 0);
                S += log_q(_em[r] + kin, _total[r] + delta);
                S += log_q(_ep[r] + kout, _total[r] + delta);
                return S;
            };

        auto get_Sr = [&](int delta)
            {
                assert(_total[r] + delta + 1 >= 0);
                return lgamma_fast(_total[r] + delta + 1);
            };

        auto get_Sk = [&](pair<size_t, size_t>& deg, int delta)
            {
                int nd = 0;
                auto iter = _hist[r].find(deg);
                if (iter != _hist[r].end())
                    nd = iter->second;
                assert(nd + delta >= 0);
                return -lgamma_fast(nd + delta + 1);
            };

692
693
694
695
696
697
698
699
700
701
702
703
704
        double S_b = 0, S_a = 0;
        int tkin = 0, tkout = 0, n = 0;
        for (auto& k : ks)
        {
            size_t kin = get<0>(k);
            size_t kout = get<1>(k);
            int nk = get<2>(k);
            tkin += kin * nk;
            if (_ignore_degree[v] != 2)
                tkout += kout * nk;
            n += nk;

            auto deg = make_pair(kin, kout);
705
706
            S_b += get_Sk(deg,         0);
            S_a += get_Sk(deg, diff * nk);
707
708
        }

709
710
        S_b += get_Se(       0,           0,            0);
        S_a += get_Se(diff * n, diff * tkin, diff * tkout);
711

712
713
        S_b += get_Sr(       0);
        S_a += get_Sr(diff * n);
714
715
716
717
718
719
720
721
722
723
724
725
726
727

        return S_a - S_b;
    }

    template <class Graph, class VWeight, class EWeight, class Degs>
    void change_vertex(size_t v, size_t r, bool deg_corr, Graph& g,
                       VWeight& vweight, EWeight& eweight, Degs& degs,
                       int diff)
    {
        auto&& ks = get_degs(v, vweight, eweight, degs, g);
        for (auto& k : ks)
        {
            auto kin = get<0>(k);
            auto kout = get<1>(k);
728
            int n = get<2>(k);
729
730
731
732
733
734
735
736
            change_k(v, r, deg_corr, n, kin, kout, diff);
        }
    }

    template <class Graph, class VWeight, class EWeight, class Degs>
    void remove_vertex(size_t v, size_t r, bool deg_corr, Graph& g,
                       VWeight& vweight, EWeight& eweight, Degs& degs)
    {
737
738
        if (r == null_group || vweight[v] == 0)
            return;
739
740
741
742
743
744
745
746
        r = get_r(r);
        change_vertex(v, r, deg_corr, g, vweight, eweight, degs, -1);
    }

    template <class Graph, class VWeight, class EWeight, class Degs>
    void add_vertex(size_t v, size_t nr, bool deg_corr, Graph& g,
                    VWeight& vweight, EWeight& eweight, Degs& degs)
    {
747
748
        if (nr == null_group || vweight[v] == 0)
            return;
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
        nr = get_r(nr);
        change_vertex(v, nr, deg_corr, g, vweight, eweight, degs, 1);
    }

    void change_k(size_t v, size_t r, bool deg_corr, int vweight,
                  int kin, int kout, int diff)
    {
        if (_total[r] == 0 && diff * vweight > 0)
            _actual_B++;

        if (_total[r] == vweight && diff * vweight < 0)
            _actual_B--;

        _total[r] += diff * vweight;
        _N += diff * vweight;

765
766
        assert(_total[r] >= 0);

767
768
769
770
771
772
773
774
775
776
777
        if (deg_corr && _ignore_degree[v] != 1)
        {
            if (_ignore_degree[v] == 2)
                kout = 0;
            auto deg = make_pair(kin, kout);
            _hist[r][deg] += diff * vweight;
            _em[r] += diff * deg.first * vweight;
            _ep[r] += diff * deg.second * vweight;
        }
    }

778
779
780
781
782
    size_t get_N()
    {
        return _N;
    }

783
784
785
786
787
    size_t get_actual_B()
    {
        return _actual_B;
    }

788
789
790
791
792
793
private:
    vector<size_t>& _bmap;
    size_t _N;
    size_t _E;
    size_t _actual_B;
    size_t _total_B;
794
    bool _allow_empty;
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
    vector<map_t> _hist;
    vector<int> _total;
    vector<int> _ep;
    vector<int> _em;
    vector<uint8_t> _ignore_degree;
};

// ===============================
// Block moves
// ===============================

// this structure speeds up the access to the edges between given blocks, since
// we're using an adjacency list to store the block structure (it is simply an
// adjacency matrix)

810
template <class BGraph>
811
812
813
class EMat
{
public:
814
815
    template <class RNG>
    EMat(BGraph& bg, RNG&)
816
    {
817
        sync(bg);
818
819
    }

820
    void sync(BGraph& bg)
821
822
823
    {
        size_t B = num_vertices(bg);
        _mat.resize(boost::extents[B][B]);
824
825
        std::fill(_mat.data(), _mat.data() + _mat.num_elements(), _null_edge);

826
827
        for (auto e : edges_range(bg))
        {
828
            assert(get_me(source(e, bg),target(e, bg)) == _null_edge);
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
            _mat[source(e, bg)][target(e, bg)] = e;
            if (!is_directed::apply<BGraph>::type::value)
                _mat[target(e, bg)][source(e, bg)] = e;
        }
    }

    typedef typename graph_traits<BGraph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<BGraph>::edge_descriptor edge_t;

    const auto& get_me(vertex_t r, vertex_t s) const
    {
        return _mat[r][s];
    }

    void put_me(vertex_t r, vertex_t s, const edge_t& e)
    {
        _mat[r][s] = e;
        if (!is_directed::apply<BGraph>::type::value && r != s)
            _mat[s][r] = e;
    }

850
    void remove_me(const edge_t& me, BGraph& bg)
851
    {
852
853
854
855
856
857
        auto r = source(me, bg);
        auto s = target(me, bg);
        _mat[r][s] = _null_edge;
        if (!is_directed::apply<BGraph>::type::value)
            _mat[s][r] = _null_edge;
        remove_edge(me, bg);
858
859
860
861
862
863
864
865
866
    }

    const auto& get_null_edge() const { return _null_edge; }

private:
    multi_array<edge_t, 2> _mat;
    static const edge_t _null_edge;
};

867
868
template <class BGraph>
const typename EMat<BGraph>::edge_t EMat<BGraph>::_null_edge;
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894


template <class Key>
class perfect_hash_t
{
public:
    template <class RNG>
    perfect_hash_t(size_t N, RNG& rng)
        : _index(std::make_shared<std::vector<size_t>>())
    {
        auto& index = *_index;
        index.reserve(N);
        for (size_t i = 0; i < N; ++i)
            index.push_back(i);
        std::shuffle(index.begin(), index.end(), rng);
    }
    perfect_hash_t() {}
    size_t operator()(const Key& k) const {return (*_index)[k];}
private:
    std::shared_ptr<std::vector<size_t>> _index;
};

// 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 above, but takes less space and is slower)

895
template <class BGraph>
896
897
898
899
class EHash
{
public:

900
901
    template <class RNG>
    EHash(BGraph& bg, RNG& rng)
902
        : _hash_function(num_vertices(bg), rng),
903
          _hash(num_vertices(bg), ehash_t(0, _hash_function))
904
    {
905
        sync(bg);
906
907
    }

908
    void sync(BGraph& bg)
909
910
911
    {
        _hash.clear();
        _hash.resize(num_vertices(bg), ehash_t(0, _hash_function));
912
913

        for (auto e : edges_range(bg))
914
915
        {
            assert(get_me(source(e, bg), target(e, bg)) == _null_edge);
916
            put_me(source(e, bg), target(e, bg), e);
917
        }
918
919
920
921
922
923
924
    }

    typedef typename graph_traits<BGraph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<BGraph>::edge_descriptor edge_t;

    const auto& get_me(vertex_t r, vertex_t s) const
    {
925
926
        if (!is_directed::apply<BGraph>::type::value && r > s)
            std::swap(r, s);
927
928
929
930
931
932
933
934
935
        auto& map = _hash[r];
        const auto& iter = map.find(s);
        if (iter == map.end())
            return _null_edge;
        return iter->second;
    }

    void put_me(vertex_t r, vertex_t s, const edge_t& e)
    {
936
937
        if (!is_directed::apply<BGraph>::type::value && r > s)
            std::swap(r, s);
938
939
940
941
        assert(r < _hash.size());
        _hash[r][s] = e;
    }

942
    void remove_me(const edge_t& me, BGraph& bg)
943
    {
944
945
946
947
948
949
950
        auto r = source(me, bg);
        auto s = target(me, bg);
        if (!is_directed::apply<BGraph>::type::value && r > s)
            std::swap(r, s);
        assert(r < _hash.size());
        _hash[r].erase(s);
        remove_edge(me, bg);
951
952
953
954
955
956
957
958
959
960
961
    }

    const auto& get_null_edge() const { return _null_edge; }

private:
    perfect_hash_t<vertex_t> _hash_function;
    typedef gt_hash_map<vertex_t, edge_t, perfect_hash_t<vertex_t>> ehash_t;
    std::vector<ehash_t> _hash;
    static const edge_t _null_edge;
};

962
963
template <class BGraph>
const typename EHash<BGraph>::edge_t EHash<BGraph>::_null_edge;
964

965
966
967
template <class Vertex, class Eprop, class Emat, class BEdge>
inline auto get_beprop(Vertex r, Vertex s, const Eprop& eprop, const Emat& emat,
                       BEdge& me)
968
{
969
970
    typedef typename property_traits<Eprop>::value_type val_t;
    me = emat.get_me(r, s);
971
    if (me != emat.get_null_edge())
972
973
974
975
976
977
978
979
980
981
        return eprop[me];
    return val_t();
}

template <class Vertex, class Eprop, class Emat>
inline auto get_beprop(Vertex r, Vertex s, const Eprop& eprop, const Emat& emat)
{
    typedef typename property_traits<Eprop>::key_type bedge_t;
    bedge_t me;
    return get_beprop(r, s, eprop, emat, me);
982
983
984
985
986
987
}


// Manage a set of block pairs and corresponding edge counts that will be
// updated

988
template <class Graph, class BGraph, class... EVals>
989
990
991
class EntrySet
{
public:
992
993
    typedef typename graph_traits<BGraph>::edge_descriptor bedge_t;

994
995
996
997
998
999
1000
    EntrySet(size_t B)
    {
        _r_field_t.resize(B, _null);
        _nr_field_t.resize(B, _null);

        if (is_directed::apply<Graph>::type::value)
        {
For faster browsing, not all history is shown. View entire blame