graph_blockmodel_util.hh 42.9 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
29
30
// 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"

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

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

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

namespace graph_tool
{

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
// 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>>{}>{});
}

80
81
82
83
// ====================
// Entropy calculation
// ====================

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
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;
};
102
103
104
105

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

106
107
108
109
110
111
112
113
114
115
116
117
// 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
    {
118
#ifndef __clang__
119
        constexpr double log_2 = log(2);
120
121
122
#else
        const double log_2 = log(2);
#endif
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
        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);
    }
}

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
// "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));
}

178
179


180
181
182
183
184
185
186
187
188
189
190
191
192
193
// 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.;

194
195
    assert(wr_r + wr_s > 0);

196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
    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;
}

216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
// Weighted entropy terms

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);
}

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);
}


242
243
244
245
// ===============
// Partition stats
// ===============

246
247
constexpr size_t null_group = std::numeric_limits<size_t>::max();

248
249
250
251
252
253
254
255
256
257
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>
258
259
260
        degs {{make_tuple(size_t(in_degreeS()(v, g, eweight)),
                          size_t(out_degreeS()(v, g, eweight)),
                          size_t(vweight[v]))}};
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
    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,
281
282
                      const Mprop& ignore_degree, std::vector<size_t>& bmap,
                      bool allow_empty)
283
        : _bmap(bmap), _N(0), _E(E), _total_B(B), _allow_empty(allow_empty)
284
285
286
    {
        for (auto v : vlist)
        {
287
288
289
            if (vweight[v] == 0)
                continue;

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
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
            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;
344
        if (_allow_empty)
345
            S += lbinom(_total_B + _N - 1, _N);
346
347
        else
            S += lbinom(_N - 1, _actual_B - 1);
348
        S += lgamma_fast(_N + 1);
349
        for (auto nr : _total)
350
            S -= lgamma_fast(nr + 1);
351
352
353
        return S;
    }

354
    double get_deg_dl_ent()
355
356
357
358
    {
        double S = 0;
        for (size_t r = 0; r < _ep.size(); ++r)
        {
359
            size_t total = 0;
360
            for (auto& k_c : _hist[r])
361
            {
362
                S -= xlogx(k_c.second);
363
364
365
                total += k_c.second;
            }
            S += xlogx(total);
366
367
368
        }
        return S;
    }
369

370
371
372
373
374
375
376
377
378
379
    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;
    }
380

381
382
383
384
385
386
387
    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]);
388

389
            size_t total = 0;
390
            for (auto& k_c : _hist[r])
391
            {
392
                S -= lgamma_fast(k_c.second + 1);
393
394
395
                total += k_c.second;
            }
            S += lgamma_fast(total + 1);
396
397
398
399
        }
        return S;
    }

400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
    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();
        }
    }

415
    template <class VProp>
416
    double get_delta_partition_dl(size_t v, size_t r, size_t nr, VProp& vweight)
417
418
419
420
    {
        if (r == nr)
            return 0;

421
422
423
424
425
426
427
        if (r != null_group)
            r = get_r(r);

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

        int n = vweight[v];
428
        if (n == 0)
429
430
431
432
433
434
        {
            if (r == null_group)
                n = 1;
            else
                return 0;
        }
435
436
437

        double S_b = 0, S_a = 0;

438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
        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);
458
459

        int dB = 0;
460
        if (r != null_group && _total[r] == n)
461
            dB--;
462
        if (nr != null_group && _total[nr] == 0)
463
464
            dB++;

465
        if ((dN != 0 || dB != 0) && !_allow_empty)
466
        {
467
468
469
470
            //S_b += lbinom_fast(std::min(_total_B, _N), _actual_B);
            S_b += lbinom_fast(_N - 1, _actual_B - 1);
            //S_a += lbinom_fast(std::min(_total_B, _N + dN), _actual_B + dB);
            S_a += lbinom_fast(_N - 1 + dN, _actual_B + dB - 1);
471
472
473
474
475
476
477
478
479
        }

        return S_a - S_b;
    }

    template <class Graph, class VProp>
    double get_delta_edges_dl(size_t v, size_t r, size_t nr, VProp& vweight,
                              Graph&)
    {
480
        if (r == nr || _allow_empty)
481
482
            return 0;

483
484
485
486
        if (r != null_group)
            r = get_r(r);
        if (nr != null_group)
            nr = get_r(nr);
487
488
489
490
491
492

        double S_b = 0, S_a = 0;

        int n = vweight[v];

        if (n == 0)
493
494
495
496
497
498
        {
            if (r == null_group)
                n = 1;
            else
                return 0;
        }
499
500

        int dB = 0;
501
        if (r != null_group && _total[r] == n)
502
            dB--;
503
        if (nr != null_group && _total[nr] == 0)
504
505
506
507
508
509
510
511
512
513
514
515
            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;
                };

516
517
            S_b += lbinom(get_x(_actual_B) + _E - 1, _E);
            S_a += lbinom(get_x(_actual_B + dB) + _E - 1, _E);
518
519
520
521
522
523
524
        }

        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,
525
                            EProp& eweight, Degs& degs, Graph& g, int kind)
526
    {
527
        if (r == nr || _ignore_degree[v] == 1 || vweight[v] == 0)
528
            return 0;
529
530
531
532
        if (r != null_group)
            r = get_r(r);
        if (nr != null_group)
            nr = get_r(nr);
533
        auto&& ks = get_degs(v, vweight, eweight, degs, g);
534
535
536
537
538
539
540
541
542
        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;
        }
543
        double dS = 0;
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
        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();
        }
567
568
569
        return dS;
    }

570
571
    template <class Ks>
    double get_delta_deg_dl_ent_change(size_t r, Ks& ks, int diff)
572
    {
573
574
575
576
577
578
579
580
581
582
        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);
            };
583

584
585
586
587
588
589
590
591
592
593
594
595
596
597
        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);
598

599
600
601
602
603
        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)
604
    {
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
        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;
    }
630
631
632


    template <class Ks>
633
    double get_delta_deg_dl_dist_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
660
661
662

        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);
            };

663
664
665
666
667
668
669
670
671
672
673
674
675
        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);
676
677
            S_b += get_Sk(deg,         0);
            S_a += get_Sk(deg, diff * nk);
678
679
        }

680
681
        S_b += get_Se(       0,           0,            0);
        S_a += get_Se(diff * n, diff * tkin, diff * tkout);
682

683
684
        S_b += get_Sr(       0);
        S_a += get_Sr(diff * n);
685
686
687
688
689
690
691
692
693
694
695
696
697
698

        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);
699
            int n = get<2>(k);
700
701
702
703
704
705
706
707
            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)
    {
708
709
        if (r == null_group || vweight[v] == 0)
            return;
710
711
712
713
714
715
716
717
        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)
    {
718
719
        if (nr == null_group || vweight[v] == 0)
            return;
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
        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;

736
737
        assert(_total[r] >= 0);

738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
        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;
        }
    }

private:
    vector<size_t>& _bmap;
    size_t _N;
    size_t _E;
    size_t _actual_B;
    size_t _total_B;
755
    bool _allow_empty;
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
    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)

771
template <class BGraph>
772
773
774
class EMat
{
public:
775
776
    template <class RNG>
    EMat(BGraph& bg, RNG&)
777
    {
778
        sync(bg);
779
780
    }

781
    void sync(BGraph& bg)
782
783
784
    {
        size_t B = num_vertices(bg);
        _mat.resize(boost::extents[B][B]);
785
786
        std::fill(_mat.data(), _mat.data() + _mat.num_elements(), _null_edge);

787
788
        for (auto e : edges_range(bg))
        {
789
            assert(get_me(source(e, bg),target(e, bg)) == _null_edge);
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
            _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;
    }

    void remove_me(vertex_t r, vertex_t s, const edge_t& me, BGraph& bg,
812
                   bool delete_edge = true)
813
814
815
    {
        if (delete_edge)
        {
816
            _mat[r][s] = _null_edge;
817
            if (!is_directed::apply<BGraph>::type::value)
818
                _mat[s][r] = _null_edge;
819
820
821
822
823
824
825
826
827
828
829
            remove_edge(me, bg);
        }
    }

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

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

830
831
template <class BGraph>
const typename EMat<BGraph>::edge_t EMat<BGraph>::_null_edge;
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857


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)

858
template <class BGraph>
859
860
861
862
class EHash
{
public:

863
864
    template <class RNG>
    EHash(BGraph& bg, RNG& rng)
865
        : _hash_function(num_vertices(bg), rng),
866
          _hash(num_vertices(bg), ehash_t(0, _hash_function))
867
    {
868
        sync(bg);
869
870
    }

871
    void sync(BGraph& bg)
872
873
874
    {
        _hash.clear();
        _hash.resize(num_vertices(bg), ehash_t(0, _hash_function));
875
876

        for (auto e : edges_range(bg))
877
878
        {
            assert(get_me(source(e, bg), target(e, bg)) == _null_edge);
879
            put_me(source(e, bg), target(e, bg), e);
880
        }
881
882
883
884
885
886
887
    }

    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
    {
888
889
        if (!is_directed::apply<BGraph>::type::value && r > s)
            std::swap(r, s);
890
891
892
893
894
895
896
897
898
        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)
    {
899
900
        if (!is_directed::apply<BGraph>::type::value && r > s)
            std::swap(r, s);
901
902
903
904
905
906
907
908
909
        assert(r < _hash.size());
        _hash[r][s] = e;
    }

    void remove_me(vertex_t r, vertex_t s, const edge_t& me, BGraph& bg,
                   bool delete_edge = true)
    {
        if (delete_edge)
        {
910
911
            if (!is_directed::apply<BGraph>::type::value && r > s)
                std::swap(r, s);
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
            assert(r < _hash.size());
            _hash[r].erase(s);
            remove_edge(me, bg);
        }
    }

    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;
};

927
928
template <class BGraph>
const typename EHash<BGraph>::edge_t EHash<BGraph>::_null_edge;
929

930
931
932
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)
933
{
934
935
    typedef typename property_traits<Eprop>::value_type val_t;
    me = emat.get_me(r, s);
936
    if (me != emat.get_null_edge())
937
938
939
940
941
942
943
944
945
946
        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);
947
948
949
950
951
952
}


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

953
template <class Graph, class BGraph, class... EVals>
954
955
956
class EntrySet
{
public:
957
958
    typedef typename graph_traits<BGraph>::edge_descriptor bedge_t;

959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
    EntrySet(size_t B)
    {
        _r_field_t.resize(B, _null);
        _nr_field_t.resize(B, _null);

        if (is_directed::apply<Graph>::type::value)
        {
            _r_field_s.resize(B, _null);
            _nr_field_s.resize(B, _null);
        }
    }

    void set_move(size_t r, size_t nr)
    {
        _rnr = make_pair(r, nr);
    }

976
    template <class... DVals>
977
    void insert_delta(size_t t, size_t s, DVals... delta)
978
    {
979
        insert_delta_imp(t, s, typename is_directed::apply<Graph>::type(),
980
                         delta...);
981
982
    }

983
    template <class... DVals>
984
    void insert_delta_imp(size_t t, size_t s, std::true_type, DVals... delta)
985
986
987
    {
        bool src = false;
        if (t != _rnr.first && t != _rnr.second)
988
        {
989
990
            std::swap(t, s);
            src = true;
991
992
        }

993
        assert(t == _rnr.first || t == _rnr.second);
994

995
996
        auto& r_field = (src) ? _r_field_s : _r_field_t;
        auto& nr_field = (src) ? _nr_field_s : _nr_field_t;
997

998
        auto& field = (_rnr.first == t) ? r_field : nr_field;
999
1000
        if (field[s] == _null)
        {