graph_blockmodel_multiflip_mcmc.hh 26.5 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-2019 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
23
24
//
// 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_MULTIFLIP_MCMC_HH
#define GRAPH_BLOCKMODEL_MULTIFLIP_MCMC_HH

#include "config.h"

#include <vector>
#include <algorithm>
25
#include <queue>
26
27

#include "graph_tool.hh"
28
#include "../support/graph_state.hh"
29
30
31
32
33
34
35
36
37
38
39
40
41
#include "graph_blockmodel_util.hh"
#include <boost/mpl/vector.hpp>

namespace graph_tool
{
using namespace boost;
using namespace std;

#define MCMC_BLOCK_STATE_params(State)                                         \
    ((__class__,&, mpl::vector<python::object>, 1))                            \
    ((state, &, State&, 0))                                                    \
    ((beta,, double, 0))                                                       \
    ((c,, double, 0))                                                          \
42
    ((d,, double, 0))                                                          \
43
    ((psingle,, double, 0))                                                    \
44
    ((psplit,, double, 0))                                                     \
45
46
47
48
    ((pmerge,, double, 0))                                                     \
    ((pmergesplit,, double, 0))                                                \
    ((nproposal, &, vector<size_t>&, 0))                                       \
    ((nacceptance, &, vector<size_t>&, 0))                                     \
49
    ((gibbs_sweeps,, size_t, 0))                                               \
50
    ((entropy_args,, entropy_args_t, 0))                                       \
51
    ((verbose,, int, 0))                                                       \
52
    ((force_move,, bool, 0))                                                   \
53
54
    ((niter,, size_t, 0))

55
enum class move_t { single = 0, split, merge, mergesplit, null };
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

template <class State>
struct MCMC
{
    GEN_STATE_BASE(MCMCBlockStateBase, MCMC_BLOCK_STATE_params(State))

    template <class... Ts>
    class MCMCBlockState
        : public MCMCBlockStateBase<Ts...>
    {
    public:
        GET_PARAMS_USING(MCMCBlockStateBase<Ts...>,
                         MCMC_BLOCK_STATE_params(State))
        GET_PARAMS_TYPEDEF(Ts, MCMC_BLOCK_STATE_params(State))

        template <class... ATs,
                  typename std::enable_if_t<sizeof...(ATs) ==
                                            sizeof...(Ts)>* = nullptr>
        MCMCBlockState(ATs&&... as)
           : MCMCBlockStateBase<Ts...>(as...),
            _g(_state._g),
77
78
79
80
            _groups(num_vertices(_state._bg)),
            _vpos(get(vertex_index_t(), _state._g),
                  num_vertices(_state._g)),
            _rpos(get(vertex_index_t(), _state._bg),
81
82
                  num_vertices(_state._bg)),
            _bnext(get(vertex_index_t(), _state._g),
83
84
                   num_vertices(_state._g)),
            _btemp(get(vertex_index_t(), _state._g),
85
                   num_vertices(_state._g))
86
87
88
89
90
91
92
        {
            _state.init_mcmc(_c,
                             (_entropy_args.partition_dl ||
                              _entropy_args.degree_dl ||
                              _entropy_args.edges_dl));
            for (auto v : vertices_range(_state._g))
            {
93
                if (_state.node_weight(v) == 0)
94
95
                    continue;
                add_element(_groups[_state._b[v]], _vpos, v);
96
                _N += _state.node_weight(v);
97
            }
98

99
            for (auto r : vertices_range(_state._bg))
100
            {
101
102
                if (_state._wr[r] == 0)
                    continue;
103
                add_element(_rlist, _rpos, r);
104
            }
105
106
107
108
109
110
111

            std::vector<move_t> moves
                = {move_t::single, move_t::split, move_t::merge,
                   move_t::mergesplit};
            std::vector<double> probs
                = {_psingle, _psplit, _pmerge, _pmergesplit};
            _move_sampler = Sampler<move_t, mpl::false_>(moves, probs);
112
113
114
115
        }

        typename state_t::g_t& _g;

116
117
118
        std::vector<std::vector<size_t>> _groups;
        typename vprop_map_t<size_t>::type::unchecked_t _vpos;
        typename vprop_map_t<size_t>::type::unchecked_t _rpos;
119
        size_t _nmoves = 0;
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
        std::vector<std::vector<std::tuple<size_t, size_t>>> _bstack;

        Sampler<move_t, mpl::false_> _move_sampler;

        void _push_b_dispatch() {}

        template <class... Vs>
        void _push_b_dispatch(const std::vector<size_t>& vs, Vs&&... vvs)
        {
            auto& back = _bstack.back();
            for (auto v : vs)
                back.emplace_back(v, _state._b[v]);
            _push_b_dispatch(std::forward<Vs>(vvs)...);
        }

        template <class... Vs>
        void push_b(Vs&&... vvs)
        {
            _bstack.emplace_back();
            _push_b_dispatch(std::forward<Vs>(vvs)...);
        }

        void pop_b()
        {
            auto& back = _bstack.back();
            for (auto& vb : back)
            {
                size_t v = get<0>(vb);
                size_t s = get<1>(vb);
                move_vertex(v, s);
            }
            _bstack.pop_back();
        }
154

155
        std::vector<size_t> _rlist;
156
        std::vector<size_t> _vs;
157

158
159
160
        typename vprop_map_t<int>::type::unchecked_t _bnext;
        typename vprop_map_t<int>::type::unchecked_t _btemp;

161
        constexpr static move_t _null_move = move_t::null;
162

163
164
        size_t _N = 0;

165
166
167
        double _dS;
        double _a;

168
        size_t node_state(size_t r)
169
        {
170
            return r;
171
172
        }

173
        constexpr bool skip_node(size_t)
174
        {
175
            return false;
176
177
        }

178
        template <bool sample_branch=true>
179
        size_t sample_new_group(size_t v, rng_t& rng)
180
        {
181
182
            _state.get_empty_block(v);
            auto t = uniform_sample(_state._empty_blocks, rng);
183

184
            auto r = _state._b[v];
185
            _state._bclabel[t] = _state._bclabel[r];
186
            if (_state._coupled_state != nullptr)
187
188
189
            {
                if constexpr (sample_branch)
                {
190
191
192
193
194
                    do
                    {
                        _state._coupled_state->sample_branch(t, r, rng);
                    }
                    while(!_state.allow_move(r, t));
195
196
197
198
199
200
201
202
203
204
                }
                else
                {
                    auto& bh = _state._coupled_state->get_b();
                    bh[t] = bh[r];
                }
                auto& hpclabel = _state._coupled_state->get_pclabel();
                hpclabel[t] = _state._pclabel[v];
            }

205
206
207
208
209
210
211
212
            if (t >= _groups.size())
            {
                _groups.resize(t + 1);
                _rpos.resize(t + 1);
            }
            assert(_state._wr[t] == 0);
            return t;
        }
213

214
215
        void move_vertex(size_t v, size_t r)
        {
216
217
218
            size_t s = _state._b[v];
            if (s == r)
                return;
219
220
221
            remove_element(_groups[s], _vpos, v);
            _state.move_vertex(v, r);
            add_element(_groups[r], _vpos, v);
222
            _nmoves++;
223
        }
224

225
226
227
228
229

        template <class RNG>
        std::tuple<double, double>
        gibbs_sweep(std::vector<size_t>& vs, size_t r, size_t s,
                    double beta, RNG& rng)
230
        {
231
232
233
234
235
236
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
264
265
266
267
268
269
270
271
272
273
274
275
276
277
            double lp = 0, dS = 0;
            std::array<double,2> p = {0,0};
            std::shuffle(vs.begin(), vs.end(), rng);
            for (auto v : vs)
            {
                size_t bv = _state._b[v];
                size_t nbv = (bv == r) ? s : r;
                double ddS;
                if (_state.virtual_remove_size(v) > 0)
                    ddS = _state.virtual_move(v, bv, nbv, _entropy_args);
                else
                    ddS = std::numeric_limits<double>::infinity();

                if (!std::isinf(beta) && !std::isinf(ddS))
                {
                    double Z = log_sum(0., -ddS * beta);
                    p[0] = -ddS * beta - Z;
                    p[1] = -Z;
                }
                else
                {
                    if (ddS < 0)
                    {
                        p[0] = 0;
                        p[1] = -std::numeric_limits<double>::infinity();
                    }
                    else
                    {
                        p[0] = -std::numeric_limits<double>::infinity();;
                        p[1] = 0;
                    }
                }

                std::bernoulli_distribution sample(exp(p[0]));
                if (sample(rng))
                {
                    move_vertex(v, nbv);
                    lp += p[0];
                    dS += ddS;
                }
                else
                {
                    lp += p[1];
                }
            }
            return {dS, lp};
        }
278

279
        template <bool forward=true, class RNG>
280
        std::tuple<double, size_t, size_t>
281
        stage_split_random(std::vector<size_t>& vs, size_t r, size_t s, RNG& rng)
282
        {
283
284
285
            std::array<size_t, 2> rt = {null_group, null_group};
            std::array<double, 2> ps;
            double dS = 0;
286
287
288
289

            std::uniform_real_distribution<> unit(0, 1);
            double p = unit(rng);

290
291
            std::shuffle(vs.begin(), vs.end(), rng);
            for (auto v : vs)
292
293
294
            {
                if (rt[0] == null_group)
                {
295
                    rt[0] = r;
296
297
                    dS += _state.virtual_move(v, _state._b[v], rt[0],
                                              _entropy_args);
298
                    move_vertex(v, rt[0]);
299
300
                    continue;
                }
301

302
303
                if (rt[1] == null_group)
                {
304
305
                    if constexpr (forward)
                        rt[1] = (s == null_group) ? sample_new_group(v, rng) : s;
306
307
                    else
                        rt[1] = s;
308
309
                    dS += _state.virtual_move(v, _state._b[v], rt[1],
                                              _entropy_args);
310
                    move_vertex(v, rt[1]);
311
312
                    continue;
                }
313

314
315
                ps[0] = log(p);
                ps[1] = log1p(-p);
316

317
                double Z = log_sum(ps[0], ps[1]);
318
                double p0 = ps[0] - Z;
319
320
321
                std::bernoulli_distribution sample(exp(p0));
                if (sample(rng))
                {
322
323
324
                    dS += _state.virtual_move(v, _state._b[v], rt[0],
                                              _entropy_args);
                    move_vertex(v, rt[0]);
325
326
327
                }
                else
                {
328
                    dS += _state.virtual_move(v, _state._b[v], rt[1],
329
                                              _entropy_args);
330
                    move_vertex(v, rt[1]);
331
332
                }
            }
333
334
            return {dS, rt[0], rt[1]};
        }
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
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
        template <bool forward=true, class RNG>
        std::tuple<double, size_t, size_t>
        stage_split_scatter(std::vector<size_t>& vs, size_t r, size_t s, RNG& rng)
        {
            std::array<size_t, 2> rt = {null_group, null_group};
            std::array<double, 2> ps;
            double dS = 0;

            if (s != null_group && _groups[s].empty())
                _state.move_vertex(_groups[r].front(), s);

            size_t t;
            if (_rlist.size() < (forward ? _N - 1 : _N))
                t = sample_new_group<false>(_groups[r].front(), rng);
            else
                t = r;

            if (s != null_group && _groups[s].empty())
                _state.move_vertex(_groups[r].front(), r);

            for (auto v : _groups[r])
            {
                dS += _state.virtual_move(v, _state._b[v], t,
                                          _entropy_args);
                move_vertex(v, t);
            }

            if constexpr (!forward)
            {
                for (auto v : _groups[s])
                {
                    dS += _state.virtual_move(v, _state._b[v], t,
                                              _entropy_args);
                    move_vertex(v, t);
                }
            }

            std::shuffle(vs.begin(), vs.end(), rng);
            for (auto v : vs)
            {
                if (rt[0] == null_group)
                {
                    rt[0] = r;
                    dS += _state.virtual_move(v, _state._b[v], rt[0],
                                              _entropy_args);
                    move_vertex(v, rt[0]);
                    continue;
                }

                if (rt[1] == null_group)
                {
                    if constexpr (forward)
                        rt[1] = (s == null_group) ? sample_new_group(v, rng) : s;
                    else
                        rt[1] = s;
                    dS += _state.virtual_move(v, _state._b[v], rt[1],
                                              _entropy_args);
                    move_vertex(v, rt[1]);
                    continue;
                }

                ps[0] = _state.virtual_move(v, _state._b[v], rt[0],
                                            _entropy_args);
                ps[1] = _state.virtual_move(v, _state._b[v], rt[1],
                                            _entropy_args);;

                double Z = log_sum(ps[0], ps[1]);
                double p0 = ps[0] - Z;
                std::bernoulli_distribution sample(exp(p0));
                if (sample(rng))
                {
                    dS += ps[0];
                    move_vertex(v, rt[0]);
                }
                else
                {
                    dS += ps[1];
                    move_vertex(v, rt[1]);
                }
            }
            return {dS, rt[0], rt[1]};
        }
418

419
420
421
422
423
        template <class RNG, bool forward=true>
        std::tuple<size_t, double, double> split(size_t r, size_t s,
                                                 RNG& rng)
        {
            auto vs = _groups[r];
424

425
426
427
428
429
            if constexpr (!forward)
                vs.insert(vs.end(), _groups[s].begin(), _groups[s].end());

            double dS;
            std::array<size_t, 2> rt;
430
431
432
433
434
435

            std::bernoulli_distribution stage_sample(.5);
            if (stage_sample(rng))
                std::tie(dS, rt[0], rt[1]) = stage_split_random<forward>(vs, r, s, rng);
            else
                std::tie(dS, rt[0], rt[1]) = stage_split_scatter<forward>(vs, r, s, rng);
436
437
438
439
440
441
442

            for (size_t i = 0; i < _gibbs_sweeps - 1; ++i)
            {
                auto ret = gibbs_sweep(vs, rt[0], rt[1],
                                       (i < _gibbs_sweeps / 2) ? 1 : _beta,
                                       rng);
                dS += get<0>(ret);
443
            }
444

445
446
447
448
449
450
451
452
453
            double lp = 0;
            if constexpr (forward)
            {
                auto ret = gibbs_sweep(vs, rt[0], rt[1], _beta, rng);
                dS += get<0>(ret);
                lp = get<1>(ret);
            }

            return {rt[1], dS, lp};
454
455
456
        }

        template <class RNG>
457
        double split_prob(size_t r, size_t s, RNG& rng)
458
        {
459
460
            auto vs = _groups[r];
            vs.insert(vs.end(), _groups[s].begin(), _groups[s].end());
461

462
463
464
465
466
467
468
469
            push_b(vs);

            for (auto v : vs)
                _btemp[v] = _state._b[v];

            split<RNG, false>(r, s, rng);

            std::shuffle(vs.begin(), vs.end(), rng);
470

471
472
            double lp1 = 0, lp2 = 0;
            for (bool swap : std::array<bool,2>({false, true}))
473
            {
474
475
476
477
                if (!swap)
                    push_b(vs);

                for (auto v : vs)
478
                {
479
480
481
482
483
484
485
                    size_t bv = _state._b[v];
                    size_t nbv = (bv == r) ? s : r;
                    double ddS;
                    if (_state.virtual_remove_size(v) > 0)
                        ddS = _state.virtual_move(v, bv, nbv, _entropy_args);
                    else
                        ddS = std::numeric_limits<double>::infinity();
486

487
488
                    if (!std::isinf(ddS))
                        ddS *= _beta;
489

490
                    double Z = log_sum(0., -ddS);
491

492
                    size_t tbv = _btemp[v];
493
                    double p;
494
                    if ((swap) ? tbv != nbv : tbv == nbv)
495
                    {
496
                        move_vertex(v, nbv);
497
                        p = -ddS - Z;
498
                    }
499
                    else
500
                    {
501
                        p = -Z;
502
                    }
503
504
505
506
507

                    if (swap)
                        lp2 += p;
                    else
                        lp1 += p;
508
509
                }

510
                if (!swap)
511
                    pop_b();
512
            }
513

514
            pop_b();
515

516
            return log_sum(lp1, lp2) - log(2);
517
518
        }

519
520
        bool allow_merge(size_t r, size_t s)
        {
521
            return _state.allow_move(r, s);
522
523
        }

524
        double merge(size_t r, size_t s)
525
526
527
        {
            double dS = 0;

528
            auto vs = _groups[r];
529

530
            for (auto v : vs)
531
            {
532
533
534
                size_t bv = _state._b[v];
                dS +=_state.virtual_move(v, bv, s, _entropy_args);
                move_vertex(v, s);
535
536
            }

537
538
539
540
541
542
543
544
545
546
547
548
549
            return dS;
        }

        template <class RNG>
        size_t sample_move(size_t r, RNG& rng)
        {
            auto s = r;
            while (s == r)
            {
                size_t v = uniform_sample(_groups[r], rng);
                s  = _state.sample_block(v, _c, 0, rng);
            }
            return s;
550
        }
551

552
553
554
555
556
557
558
559
560
561
562
        double get_move_prob(size_t r, size_t s)
        {
            double prs = 0, prr = 0;
            for (auto v : _groups[r])
            {
                prs += _state.get_move_prob(v, r, s, _c, 0, false);
                prr += _state.get_move_prob(v, r, r, _c, 0, false);
            }
            prs /= _groups[r].size();
            prr /= _groups[r].size();
            return prs/(1-prr);
563
564
        }

565
        double merge_prob(size_t r, size_t s)
566
        {
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
            return log(get_move_prob(r, s));
        }

        template <class RNG>
        std::tuple<size_t, double, double, double>
        sample_merge(size_t r, RNG& rng)
        {
            size_t s = sample_move(r, rng);

            if (!allow_merge(r, s))
                return {null_group, 0., 0., 0.};

            double pf = 0, pb = 0;
            if (!std::isinf(_beta))
            {
                pf = merge_prob(r, s);
                pb = split_prob(s, r, rng);
            }

            if (_verbose)
                cout << "merge " << _groups[r].size() << " " << _groups[s].size();

            double dS = merge(r, s);

            if (_verbose)
                cout << " " << dS << " " << pf << "  " << pb << endl;

            return {s, dS, pf, pb};
        }

        template <class RNG>
        std::tuple<size_t, double, double, double>
        sample_split(size_t r, size_t s, RNG& rng)
        {
            double dS, pf, pb=0;
            std::tie(s, dS, pf) = split(r, s, rng);
            if (!std::isinf(_beta))
                pb = merge_prob(s, r);

            if (_verbose)
                cout << "split " << _groups[r].size() << " " << _groups[s].size()
                     << " " << dS << " " << pf << " " << pb << endl;

            return {s, dS, pf, pb};
611
612
613
        }

        template <class RNG>
614
        std::tuple<move_t,size_t> move_proposal(size_t r, RNG& rng)
615
616
        {
            double pf = 0, pb = 0;
617
            _dS = _a = 0;
618
619
            _vs.clear();
            _nmoves = 0;
620

621
622
623
            move_t move = _move_sampler.sample(rng);

            switch (move)
624
            {
625
            case move_t::single:
626
                {
627
628
629
                    auto v = uniform_sample(_groups[r], rng);
                    auto s = _state.sample_block(v, _c, _d, rng);
                    if (s >= _groups.size())
630
                    {
631
632
                        _groups.resize(s+ 1);
                        _rpos.resize(s+ 1);
633
                    }
634
635
636
637
638
639
                    if (r == s || !_state.allow_move(r, s))
                        return {_null_move, 1};
                    if (_d == 0 && _groups[r].size() == 1 && !std::isinf(_beta))
                        return {_null_move, 1};
                    _dS = _state.virtual_move(v, r, s, _entropy_args);
                    if (!std::isinf(_beta))
640
                    {
641
642
643
644
645
646
647
648
649
650
651
652
                        pf = log(_state.get_move_prob(v, r, s, _c, _d, false));
                        pb = log(_state.get_move_prob(v, s, r, _c, _d, true));

                        pf += -safelog_fast(_rlist.size());
                        pf += -safelog_fast(_groups[r].size());
                        int dB = 0;
                        if (_groups[s].empty())
                            dB++;
                        if (_groups[r].size() == 1)
                            dB--;
                        pb += -safelog_fast(_rlist.size() + dB);
                        pb += -safelog_fast(_groups[s].size() + 1);
653
                    }
654
655
656
657
                    _vs.clear();
                    _vs.push_back(v);
                    _bnext[v] = s;
                    _nmoves++;
658
                }
659
660
661
                break;

            case move_t::split:
662
                {
663
                    if (_groups[r].size() < 2)
664
                        return {_null_move, 1};
665
666
667
668
669
670
671
672

                    _vs = _groups[r];
                    push_b(_vs);

                    size_t s;
                    std::tie(s, _dS, pf) = split(r, null_group, rng);

                    if (!std::isinf(_beta))
673
                    {
674
675
                        pf += log(_psplit);
                        pb = merge_prob(s, r) + log(_pmerge);
676
                    }
677

678
679
680
681
                    if (_verbose)
                        cout << "split proposal: " << _groups[r].size() << " "
                             << _groups[s].size() << " " << _dS << " " << pb - pf
                             << " " << -_dS + pb - pf << endl;
682

683
684
685
686
687
                    for (auto v : _vs)
                        _bnext[v] = _state._b[v];
                    pop_b();
                }
                break;
688

689
690
691
692
            case move_t::merge:
                {
                    if (_rlist.size() == 1)
                        return {_null_move, 1};
693

694
695
696
697
                    auto s = sample_move(r, rng);

                    if (!allow_merge(r, s))
                        return {_null_move, 1};
698
699
700

                    if (!std::isinf(_beta))
                    {
701
702
                        pf = merge_prob(r, s) + log(_pmerge);
                        pb = split_prob(s, r, rng) + log(_psplit);
703
704
                    }

705
706
707
708
709
710
711
712
713
                    _vs = _groups[r];
                    push_b(_vs);

                    _dS = merge(r, s);

                    for (auto v : _vs)
                        _bnext[v] = _state._b[v];
                    pop_b();

714
                    if (_verbose)
715
716
                        cout << "merge proposal: " <<  _groups[r].size() << " "
                             << _groups[s].size() << " " << _dS << " " << pb - pf
717
                             << " " << -_dS + pb - pf << endl;
718
                }
719
720
721
                break;

            case move_t::mergesplit:
722
                {
723
724
725
726
727
728
729
730
731
                    if (_rlist.size() == 1)
                        return {_null_move, 1};

                    push_b(_groups[r]);

                    auto ret = sample_merge(r, rng);
                    size_t s = get<0>(ret);

                    if (s == null_group)
732
733
734
                    {
                        while (!_bstack.empty())
                            pop_b();
735
                        return {_null_move, 1};
736
                    }
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762

                    _dS += get<1>(ret);
                    pf += get<2>(ret);
                    pb += get<3>(ret);

                    push_b(_groups[s]);

                    ret = sample_split(s, r, rng);
                    _dS += get<1>(ret);
                    pf += get<2>(ret);
                    pb += get<3>(ret);

                    for (auto& vs : _bstack)
                        for (auto& vb : vs)
                        {
                            auto v = get<0>(vb);
                            _vs.push_back(v);
                            _bnext[v] = _state._b[v];
                        }

                    while (!_bstack.empty())
                        pop_b();

                    if (_verbose)
                        cout << "mergesplit proposal: " << _dS << " " << pb - pf
                             << " " << -_dS + pb - pf << endl;
763
                }
764
765
766
767
768
                break;

            default:
                return {_null_move, 0};
                break;
769
            }
770

771
            _a = pb - pf;
772

773
774
775
776
777
778
            if (size_t(move) >= _nproposal.size())
            {
                _nproposal.resize(size_t(move) + 1);
                _nacceptance.resize(size_t(move) + 1);
            }
            _nproposal[size_t(move)]++;
779

780
781
782
783
784
785
786
            if (_force_move)
            {
                _nmoves = std::numeric_limits<size_t>::max();
                _a = _dS * _beta + 1;
            }

            return {move, _nmoves};
787
788
789
790
791
792
        }

        std::tuple<double, double>
        virtual_move_dS(size_t, move_t)
        {
            return {_dS, _a};
793
794
        }

795
        void perform_move(size_t, move_t move)
796
        {
797
            for (auto v : _vs)
798
            {
799
800
801
802
                size_t r = _state._b[v];
                size_t s = _bnext[v];
                if (_groups[s].empty())
                    add_element(_rlist, _rpos, s);
803

804
                move_vertex(v, s);
805

806
                if (_groups[r].empty() && has_element(_rlist, _rpos, r))
807
                    remove_element(_rlist, _rpos, r);
808
            }
809
810

            _nacceptance[size_t(move)]++;
811
        }
812

813
        constexpr bool is_deterministic()
814
        {
815
            return false;
816
817
        }

818
        constexpr bool is_sequential()
819
        {
820
            return false;
821
822
823
824
        }

        auto& get_vlist()
        {
825
826
827
828
829
830
            return _rlist;
        }

        size_t get_N()
        {
            return _N;
831
832
833
834
835
836
837
838
839
        }

        double get_beta()
        {
            return _beta;
        }

        size_t get_niter()
        {
840
            return _niter;
841
842
        }

843
        constexpr void step(size_t, move_t)
844
845
        {
        }
846
847
848
    };
};

849
std::ostream& operator<<(std::ostream& os, move_t move);
850
851
852
853

} // graph_tool namespace

#endif //GRAPH_BLOCKMODEL_MCMC_HH