graph_blockmodel_layers.hh 40 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-2018 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//
// 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_LAYERS_HH
#define GRAPH_BLOCKMODEL_LAYERS_HH

21
22
#define GRAPH_BLOCKMODEL_LAYERS_ENABLE

23
24
25
26
#include "config.h"

#include <vector>

27
28
#include "../support/graph_state.hh"
#include "../blockmodel/graph_blockmodel_util.hh"
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include "graph_blockmodel_layers_util.hh"

namespace graph_tool
{
using namespace boost;
using namespace std;

typedef eprop_map_t<int32_t>::type emap_t;
typedef vprop_map_t<std::vector<int32_t>>::type vcvmap_t;

typedef gt_hash_map<size_t, size_t> bmap_t;
typedef std::vector<bmap_t> vbmap_t;

#define LAYERED_BLOCK_STATE_params                                             \
    ((__class__,&, mpl::vector<python::object>, 1))                            \
    ((layer_states,, python::object, 0))                                       \
    ((ec,, emap_t, 0))                                                         \
    ((vc,, vcvmap_t, 0))                                                       \
    ((vmap,, vcvmap_t, 0))                                                     \
    ((block_map, &, vbmap_t&, 0))                                              \
    ((master,, bool, 0))

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

class LayeredBlockStateVirtualBase
    : public BlockStateVirtualBase
{
public:
    virtual BlockStateVirtualBase& get_layer(size_t l) = 0;
    virtual size_t get_block(size_t l, size_t v) = 0;
    virtual void set_block(size_t l, size_t v, size_t r) = 0;
    virtual size_t get_vweight(size_t l, size_t v) = 0;
    virtual void add_layer_node(size_t l, size_t v, size_t u) = 0;
    virtual void remove_layer_node(size_t l, size_t v, size_t u) = 0;
    virtual size_t get_layer_node(size_t l, size_t v) = 0;
    virtual size_t get_block_map(size_t l, size_t r, bool put_new=true) = 0;
    virtual bool check_layers() = 0;
};

67
68
69
70
71
72
73
74
template <class BaseState>
struct Layers
{
    GEN_STATE_BASE(LayeredBlockStateBase, LAYERED_BLOCK_STATE_params)

    template <class... Ts>
    class LayeredBlockState
        : public LayeredBlockStateBase<Ts...>,
75
76
          public BaseState,
          public LayeredBlockStateVirtualBase
77
78
79
80
81
82
83
84
85
86
87
88
    {
    public:
        GET_PARAMS_USING(LayeredBlockStateBase<Ts...>,
                         LAYERED_BLOCK_STATE_params)
        GET_PARAMS_TYPEDEF(Ts, LAYERED_BLOCK_STATE_params)

        GET_PARAMS_USING(BaseState, BASE_STATE_params)
        using BaseState::_bg;
        using BaseState::_m_entries;
        using BaseState::_emat;
        using BaseState::_partition_stats;
        using BaseState::is_partition_stats_enabled;
89
        using BaseState::get_move_entries;
90
91
92
93
94
95
96

        typedef vprop_map_t<int32_t>::type block_rmap_t;

        class LayerState
            : public BaseState
        {
        public:
97
98
99
            LayerState(const BaseState& base_state, LayeredBlockState& lstate,
                       bmap_t& block_map, block_rmap_t block_rmap,
                       vector<size_t>& free_blocks, size_t l)
100
                : BaseState(base_state),
101
                  _lstate(&lstate),
102
103
104
                  _block_map(block_map),
                  _block_rmap(block_rmap),
                  _free_blocks(free_blocks),
105
                  _l(l), _E(0)
106
107
108
109
110
            {
                for (auto e : edges_range(BaseState::_g))
                    _E += BaseState::_eweight[e];
            }

111
            LayeredBlockState* _lstate;
112
113
114
            bmap_t& _block_map;
            block_rmap_t _block_rmap;
            vector<size_t>& _free_blocks;
115
            size_t _l;
116
117
            size_t _E;

118
119
120
121
            using BaseState::_bg;
            using BaseState::_wr;
            using BaseState::add_block;

122
123
124
125
126
127
128
            size_t get_block_map(size_t r, bool put_new=true)
            {
                size_t r_u;
                auto iter = _block_map.find(r);
                if (iter == _block_map.end())
                {
                    if (_free_blocks.empty())
129
130
                        _free_blocks.push_back(_block_map.size());
                    r_u = _free_blocks.back();
131
                    while (r_u >= num_vertices(_bg))
132
                        add_block();
133
                    assert(r_u < num_vertices(_bg));
134

135
136
137
138
                    if (put_new)
                    {
                        _block_map[r] = r_u;
                        _block_rmap[r_u] = r;
139
140
141
                        if (_lstate->_lcoupled_state != nullptr)
                        {
                            _lstate->_lcoupled_state->add_layer_node(_l, r, r_u);
142
143
144
145
146
147
148
                            auto& hb = _lstate->_lcoupled_state->get_b();
                            auto& hb_u = BaseState::_coupled_state->get_b();
                            hb_u[r_u] = _lstate->_lcoupled_state->get_block_map(_l, hb[r]);
                            // assert(_lstate->_lcoupled_state->get_vweight(_l, r_u) == (_wr[r_u] > 0));
                            // if (_wr[r_u] == 0)
                            //     _lstate->_lcoupled_state->set_block(_l, r_u, hb[r_u]);
                            // assert(_lstate->_lcoupled_state->get_block(_l, r_u) == size_t(hb[r_u]));
149
150
151
152
                        }
                        _free_blocks.pop_back();
                        assert(_lstate->_lcoupled_state == nullptr ||
                               r_u == _lstate->_lcoupled_state->get_layer_node(_l, r));
153
154
155
156
157
158
159
160
161
162
163
164
165
                        // assert(_lstate->_lcoupled_state == nullptr ||
                        //        size_t(_bclabel[r_u]) ==
                        //        _lstate->_lcoupled_state->
                        //        get_block_map(_l, _lstate->_bclabel[r], false));
                    }
                    else
                    {
                        if (_lstate->_lcoupled_state != nullptr)
                        {
                            auto& hb = _lstate->_lcoupled_state->get_b();
                            auto& hb_u = BaseState::_coupled_state->get_b();
                            hb_u[r_u] = _lstate->_lcoupled_state->get_block_map(_l, hb[r], false);
                        }
166
167
168
169
170
                    }
                }
                else
                {
                    r_u = iter->second;
171
172
173
                    assert(size_t(_block_rmap[r_u]) == r);
                    assert(_lstate->_lcoupled_state == nullptr ||
                           r_u == _lstate->_lcoupled_state->get_layer_node(_l, r));
174
                }
175
176
177
178
179
180
181
182

                if (_lstate->_lcoupled_state != nullptr)
                {
                    auto& hb = _lstate->_lcoupled_state->get_b();
                    auto& hb_u = BaseState::_coupled_state->get_b();
                    hb_u[r_u] = _lstate->_lcoupled_state->get_block_map(_l, hb[r], put_new);
                }

183
                assert(r_u < num_vertices(_bg));
184
185
186
187
188
189
190
191
192
193
194
195
196
                return r_u;
            }

            bool has_block_map(size_t r)
            {
                return _block_map.find(r) != _block_map.end();
            }
        };

        template <class... ATs,
                  typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
        LayeredBlockState(const BaseState& base_state, ATs&&... args)
            : LayeredBlockStateBase<Ts...>(std::forward<ATs>(args)...),
197
198
              BaseState(base_state), _vc_c(_vc.get_checked()),
              _vmap_c(_vmap.get_checked())
199
        {
200
            for (int l = 0; l < python::len(_layer_states); ++l)
201
            {
202
                auto ostate = _layer_states[l];
203
204
205
206
207
208
                BaseState& state = python::extract<BaseState&>(ostate.attr("_state"));
                boost::python::object temp = ostate.attr("block_rmap").attr("_get_any")();
                boost::any& a = python::extract<boost::any&>(temp);
                block_rmap_t block_rmap = boost::any_cast<block_rmap_t>(a);
                std::vector<size_t>& free_blocks =
                    python::extract<std::vector<size_t>&>(ostate.attr("free_blocks"));
209
210
                bmap_t& block_map = _block_map[l];
                _layers.emplace_back(state, *this, block_map, block_rmap, free_blocks, l);
211
212
            }
            for (auto r : vertices_range(BaseState::_bg))
213
            {
214
                if (BaseState::_wr[r] > 0)
215
                    _actual_B++;
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
            }
            _N = BaseState::get_N();
            // assert(check_layers());
            // assert(check_edge_counts());
        }

        LayeredBlockState(const LayeredBlockState& other)
            : LayeredBlockStateBase<Ts...>(static_cast<const LayeredBlockStateBase<Ts...>&>(other)),
              BaseState(other),
              _layers(other._layers),
              _actual_B(other._actual_B),
              _N(other._N),
              _is_partition_stats_enabled(other._is_partition_stats_enabled),
              _lcoupled_state(other._lcoupled_state),
              _vc_c(_vc.get_checked()),
            _vmap_c(_vmap.get_checked())
        {
            for (auto& state : _layers)
                state._lstate = this;
235
236
237
        }

        std::vector<LayerState> _layers;
238
239
240
241
242
243
244
        size_t _actual_B = 0;
        size_t _N = 0;
        bool _is_partition_stats_enabled = false;
        LayeredBlockStateVirtualBase* _lcoupled_state = nullptr;
        typename vc_t::checked_t _vc_c;
        typename vmap_t::checked_t _vmap_c;
        openmp_mutex _llock;
245
246
247

        void move_vertex(size_t v, size_t s)
        {
248
            // assert(check_layers());
249
            // assert(check_edge_counts());
250

251
252
253
254
255
256
257
258
259
260
261
262
            if (BaseState::_vweight[v] == 0)
            {
                _b[v] = s;
                return;
            }

            size_t r = _b[v];

            if (r == s)
                return;

            if (_wr[s] == 0)
263
                _bclabel[s] = _bclabel[r];
264

265
266
            assert((_bclabel[r] == _bclabel[s]));

267
268
269
270
271
272
273
274
275
            auto& ls = _vc[v];
            auto& vs = _vmap[v];
            for (size_t j = 0; j < ls.size(); ++j)
            {
                int l = ls[j];
                size_t u = vs[j];

                auto& state = _layers[l];

276
277
                if (state._vweight[u] == 0)
                    continue;
278

279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
                assert(state.has_block_map(r));
                assert(size_t(state._b[u]) == state.get_block_map(r, false));
                assert(_lcoupled_state == nullptr ||
                       _lcoupled_state->get_vweight(l, state._b[u]) > 0);
                assert(state._wr[state._b[u]] > 0);
                size_t s_u = state.get_block_map(s);

                assert(size_t(state._b[u]) != s_u);

                state.move_vertex(u, s_u);

                assert(state._wr[s_u] > 0);
                assert(s_u == state.get_block_map(s, false));
            }

            // bottom update needs to be last, due to _coupled_state, and the
            // fact that the upper levels are affected by get_block_map()

            if (_wr[s] == 0)
                _actual_B++;

300
            // BaseState::check_edge_counts();
301
            BaseState::move_vertex(v, s);
302
            // BaseState::check_edge_counts();
303
304
305
306
307
308
309

            if (_wr[r] == 0)
                _actual_B--;

            if (_lcoupled_state != nullptr)
            {
                for (size_t j = 0; j < ls.size(); ++j)
310
                {
311
312
313
314
315
316
317
318
319
320
321
322
323
324
                    int l = ls[j];
                    size_t u = vs[j];
                    auto& state = _layers[l];

                    if (state._vweight[u] == 0)
                        continue;

                    size_t r_u = state._b[u];
                    assert(r_u == state.get_block_map(s));
                    assert(state._wr[r_u] > 0);

                    _lcoupled_state->get_layer(l).set_vertex_weight(r_u, 1);

                    r_u = state.get_block_map(r);
325
                    if (state._wr[r_u] == 0)
326
327
                        _lcoupled_state->get_layer(l).set_vertex_weight(r_u, 0);
                    assert(state._wr[r_u] == 0 || BaseState::_wr[r] != 0);
328
329
                }
            }
330
331

            // assert(check_layers());
332
            // assert(check_edge_counts());
333
334
        }

335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
        template <class Vec>
        void move_vertices(Vec& v, Vec& nr)
        {
            for (size_t i = 0; i < std::min(v.size(), nr.size()); ++i)
                move_vertex(v[i], nr[i]);
        }

        void move_vertices(python::object ovs, python::object ors)
        {
            multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
            multi_array_ref<uint64_t, 1> rs = get_array<uint64_t, 1>(ors);
            if (vs.size() != rs.size())
                throw ValueException("vertex and group lists do not have the same size");
            move_vertices(vs, rs);
        }

351
352
353
354
355
356
357
358
359
360
361
362
363
        void remove_vertex(size_t v)
        {
            size_t r = _b[v];
            auto& ls = _vc[v];
            auto& vs = _vmap[v];
            for (size_t j = 0; j < ls.size(); ++j)
            {
                int l = ls[j];
                size_t u = vs[j];
                auto& state = _layers[l];
                state.remove_vertex(u);
            }
            BaseState::remove_vertex(v);
364
365
            if (_wr[r] == 0)
                _actual_B--;
366
367
        }

368
369
370
371
        template <class Vec>
        void remove_vertices(Vec& vs)
        {
            gt_hash_map<size_t, vector<size_t>> lvs;
372
            gt_hash_set<size_t> rs;
373
            for (auto v : vs)
374
            {
375
376
                for (auto l : _vc[v])
                    lvs[l].push_back(v);
377
378
                rs.insert(_b[v]);
            }
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
            for (auto& lv : lvs)
            {
                auto l = lv.first;
                auto& state = _layers[l];
                vector<size_t> us;
                gt_hash_map<size_t, size_t> rus;
                for (auto v : lv.second)
                {
                    auto u = _vmap[v][l];
                    us.push_back(u);
                    size_t r = _b[v];
                    size_t r_u = state._b[u];
                    rus[r] = r_u;
                }
                state.remove_vertices(us);

395
396
397
398
399
                // for (auto rr_u : rus)
                // {
                //     if (state._wr[rr_u.second] == 0)
                //         state.remove_block_map(rr_u.first);
                // }
400
401
            }
            BaseState::remove_vertices(vs);
402
403
404
405
406
            for (auto r : rs)
            {
                if (_wr[r] == 0)
                    _actual_B--;
            }
407
408
409
410
411
412
413
414
        }

        void remove_vertices(python::object ovs)
        {
            multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
            remove_vertices(vs);
        }

415
416
417
418
419
420
421
422
423
424
425
426
        void add_vertex(size_t v, size_t r)
        {
            auto& ls = _vc[v];
            auto& vs = _vmap[v];
            for (size_t j = 0; j < ls.size(); ++j)
            {
                int l = ls[j];
                size_t u = vs[j];
                auto& state = _layers[l];
                size_t r_u = state.get_block_map(r);
                state.add_vertex(u, r_u);
            }
427
428
            if (_wr[r] == 0)
                _actual_B++;
429
430
431
            BaseState::add_vertex(v, r);
        }

432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
        template <class Vs, class Rs>
        void add_vertices(Vs& vs, Rs& rs)
        {
            if (vs.size() != rs.size())
                throw ValueException("vertex and group lists do not have the same size");

            gt_hash_map<size_t, vector<size_t>> lvs;
            gt_hash_map<size_t, size_t> vrs;
            for (size_t i = 0; i < vs.size(); ++i)
            {
                auto v = vs[i];
                vrs[v] = rs[i];
                for (auto l : _vc[v])
                    lvs[l].push_back(v);
            }

            for (auto& lv : lvs)
            {
                auto l = lv.first;
                auto& state = _layers[l];
                vector<size_t> us;
                vector<size_t> rus;
                for (auto v : lv.second)
                {
                    us.emplace_back(_vmap[v][l]);
                    rus.emplace_back(state.get_block_map(vrs[v]));
                }
                state.add_vertices(us, rus);
            }
461
462
463
464
465
            for (auto r : rs)
            {
                if (_wr[r] == 0)
                    _actual_B++;
            }
466
467
468
469
470
471
472
473
474
475
            BaseState::add_vertices(vs, rs);
        }

        void add_vertices(python::object ovs, python::object ors)
        {
            multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
            multi_array_ref<uint64_t, 1> rs = get_array<uint64_t, 1>(ors);
            add_vertices(vs, rs);
        }

476
477
478
479
480
481
482
483
484
485
486
487
488
489
        template <class VMap>
        void set_partition(VMap&& b)
        {
            for (auto v : vertices_range(_g))
                LayeredBlockState::move_vertex(v, b[v]);
        }

        void set_partition(boost::any& ab)
        {
            typename BaseState::b_t::checked_t& b
                = boost::any_cast<typename BaseState::b_t::checked_t&>(ab);
            set_partition(b.get_unchecked());
        }

490
        bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true)
491
        {
492
            return BaseState::allow_move(v, r, nr, allow_empty);
493
494
        }

495
        template <class MEntries>
496
        double virtual_move(size_t v, size_t r, size_t s, entropy_args_t ea,
497
498
499
                            MEntries& m_entries)
        {
            if (s == r)
500
501
            {
                m_entries.set_move(r, s, num_vertices(BaseState::_bg));
502
                return 0;
503
504
            }

505
            // assert(check_layers());
506
507
508

            double dS = 0;

509
510
511
512
513
            entropy_args_t mea(ea);
            mea.edges_dl = false;
            mea.recs = false;

            if (!_master)
514
            {
515
516
                mea.adjacency = false;
                mea.degree_dl = false;
517
            }
518
519
520
521

            dS += BaseState::virtual_move(v, r, s, mea, m_entries);

            if (_master)
522
            {
523
524
                if (ea.adjacency)
                    dS -= virtual_move_covariate(v, r, s, *this, m_entries, false);
525
526

                if (ea.edges_dl)
527
                    dS += ea.beta_dl * get_delta_edges_dl(v, r, s);
528
529
            }

530
            // assert(check_layers());
531

532
            if (ea.adjacency || ea.recs || ea.edges_dl || _lcoupled_state != nullptr)
533
            {
534
535
                scoped_lock lck(_llock);

536
537
                entropy_args_t lea(ea);
                lea.partition_dl = false;
538

539
540
541
542
543
544
545
                if (_master)
                {
                    lea.adjacency = false;
                    lea.degree_dl = false;
                    lea.edges_dl = false;
                }

546
547
548
549
550
551
                auto& ls = _vc[v];
                auto& vs = _vmap[v];
                for (size_t j = 0; j < ls.size(); ++j)
                {
                    size_t l = ls[j];
                    size_t u = vs[j];
552

553
                    auto& state = _layers[l];
554

555
556
557
                    if (state._vweight[u] == 0)
                        continue;

558
559
560
561
562
                    size_t s_u = (s != null_group) ?
                        state.get_block_map(s, false) : null_group;
                    size_t r_u = (r != null_group) ?
                        state._b[u] : null_group;

563
564
565
566
                    assert(r == null_group || state.has_block_map(r));
                    assert(r == null_group || r_u == state.get_block_map(r, false));

                    if (_master && ea.adjacency)
567
568
                        dS += virtual_move_covariate(u, r_u, s_u, state,
                                                     m_entries, true);
569

570
                    dS += state.virtual_move(u, r_u, s_u, lea, m_entries);
571
                }
572
            }
573
574
575

            // assert(check_layers());

576
577
578
            return dS;
        }

579
        double virtual_move(size_t v, size_t r, size_t s, entropy_args_t ea)
580
        {
581
            return virtual_move(v, r, s, ea, _m_entries);
582
583
        }

584
        template <class MEntries>
585
        double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
586
587
588
589
590
                             bool reverse, MEntries& m_entries)
        {
            // m_entries may include entries from different levels
            if (!reverse)
                BaseState::get_move_entries(v, r, s, m_entries);
591
            return BaseState::get_move_prob(v, r, s, c, d, reverse, m_entries);
592
593
        }

594
595
596
597
598
599
600
        double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
                             bool reverse,
                             std::vector<std::tuple<size_t, size_t, int>>& p_entries)
        {
            return BaseState::get_move_prob(v, r, s, c, d, reverse, p_entries);
        }

601
        double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
602
603
                             bool reverse)
        {
604
            return BaseState::get_move_prob(v, r, s, c, d, reverse);
605
606
        }

607
608
609
610
611
        size_t sample_block(size_t v, double c, double d, rng_t& rng)
        {
            return BaseState::sample_block(v, c, d, rng);
        }

612
613
        void merge_vertices(size_t u, size_t v)
        {
614
615
616
617
618
619
            if (u == v)
                return;

            assert(BaseState::_vweight[v] > 0);
            assert(BaseState::_vweight[u] > 0);

620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
            std::set<size_t> ls;
            gt_hash_map<size_t, size_t> ls_u, ls_v;
            for (size_t i = 0; i < _vc[u].size(); ++i)
            {
                size_t l = _vc[u][i];
                ls_u[l] = _vmap[u][i];
                ls.insert(l);
            }

            for (size_t i = 0; i < _vc[v].size(); ++i)
            {
                size_t l = _vc[v][i];
                ls_v[l] = _vmap[v][i];
                ls.insert(l);
            }

            _vc[u].clear();
            _vmap[u].clear();
            _vc[v].clear();
            _vmap[v].clear();

            for (auto l : ls)
            {
                auto iter_u = ls_u.find(l);
                auto iter_v = ls_v.find(l);

                size_t uu = (iter_u != ls_u.end()) ? iter_u->second : iter_v->second;
                size_t vv = (iter_v != ls_v.end()) ? iter_v->second : iter_u->second;

649
650
651
652
653
654
655
656
657
658
659
                auto& state = _layers[l];
                assert(state._vweight[vv] > 0);
                assert(state._vweight[uu] > 0);

                state.merge_vertices(uu, vv);

                assert(state._b[uu] == state._b[vv]);
                assert(size_t(state._b[vv]) == state.get_block_map(_b[v], false));

                assert(state._vweight[uu] > 0 || total_degreeS()(uu, state._g, state._eweight) == 0);
                assert(state._vweight[vv] > 0 || total_degreeS()(vv, state._g, state._eweight) == 0);
660
661
662
663
664

                _vc[v].push_back(l);
                _vmap[v].push_back(vv);
            }

665
666
            BaseState::merge_vertices(u, v, _ec.get_checked());

667
            assert(check_layers());
668
            // assert(check_edge_counts());
669
670
        }

671
        double entropy(entropy_args_t ea, bool propagate=false)
672
        {
673
            double S = 0, S_dl = 0;
674
675
            if (_master)
            {
676
677
678
679
680
681
682
683
                entropy_args_t mea(ea);
                mea.edges_dl = false;
                mea.recs = false;
                mea.recs_dl = false;

                S += BaseState::entropy(mea);

                if (ea.adjacency)
684
                {
685
                    S -= covariate_entropy(_bg, _mrs);
686
                    if (ea.multigraph)
687
688
689
690
                        S -= BaseState::get_parallel_entropy();
                    for (auto& state : _layers)
                    {
                        S += covariate_entropy(state._bg, state._mrs);
691
                        if (ea.multigraph)
692
693
694
695
                            S += state.get_parallel_entropy();
                    }
                }

696
                if (ea.edges_dl)
697
                {
698
699
700
                    size_t actual_B = _actual_B;
                    if (BaseState::_allow_empty)
                        actual_B = num_vertices(BaseState::_bg);
701
                    for (auto& state : _layers)
702
703
704
705
706
707
708
709
                        S_dl += get_edges_dl(actual_B, state._E, _g);
                }

                if (ea.recs)
                {
                    entropy_args_t mea = {false, false, false, false, true,
                                          false, false, false,
                                          ea.degree_dl_kind, false, ea.recs_dl,
710
                                          ea.beta_dl, false};
711
712
                    for (auto& state : _layers)
                        S += state.entropy(mea);
713
714
715
716
                }
            }
            else
            {
717
718
719
                entropy_args_t mea(ea);
                mea.partition_dl = false;
                mea.edges_dl = false;
720

721
                for (auto& state : _layers)
722
                    S += state.entropy(mea);
723

724
725
                if (ea.partition_dl)
                    S_dl += BaseState::get_partition_dl();
726

727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
                if (ea.edges_dl)
                {
                    for (auto& state : _layers)
                    {
                        size_t actual_B = 0;
                        if (BaseState::_allow_empty)
                        {
                            actual_B = num_vertices(_bg);
                        }
                        else
                        {
                            for (auto r : vertices_range(state._bg))
                                if (state._wr[r] > 0)
                                    actual_B++;
                        }
                        S_dl += get_edges_dl(actual_B, state._E, _g);
                    }
                }
745
                int L = _layers.size();
746
                S_dl += _N * (L * std::log(2) + std::log1p(-std::pow(2., -L)));
747
            }
748

749
750
751
            if (BaseState::_coupled_state != nullptr && propagate)
                S_dl += BaseState::_coupled_state->entropy(BaseState::_coupled_entropy_args, true);

752
            return S + S_dl * ea.beta_dl;
753
754
        }

755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
        double get_delta_edges_dl(size_t v, size_t r, size_t s)
        {
            if (r == s || BaseState::_allow_empty)
                return 0;
            if (BaseState::_vweight[v] == 0)
                return 0;
            int dB = 0;
            if (r != null_group && BaseState::virtual_remove_size(v) == 0)
                --dB;
            if (s != null_group && _wr[s] == 0)
                ++dB;
            double S_a = 0, S_b = 0;
            if (dB != 0)
            {
                auto get_x = [](size_t B)
                    {
771
                        if (is_directed_::apply<typename BaseState::g_t>::type::value)
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
                            return B * B;
                        else
                            return (B * (B + 1)) / 2;
                    };

                for (auto& state : _layers)
                {
                    S_b += lbinom(get_x(_actual_B) + state._E - 1, state._E);
                    S_a += lbinom(get_x(_actual_B + dB) + state._E - 1, state._E);
                }
            }
            return S_a - S_b;
        }

        double get_deg_dl(int kind)
787
788
789
        {
            if (_master)
            {
790
                return BaseState::get_deg_dl(kind);
791
792
793
794
795
            }
            else
            {
                double S = 0;
                for (auto& state : _layers)
796
                    S += state.get_deg_dl(kind);
797
798
799
800
                return S;
            }
        }

801
        double edge_entropy_term(size_t, size_t, entropy_args_t) { return 0; }
802

803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
        void enable_partition_stats()
        {
            if (!_is_partition_stats_enabled)
            {
                BaseState::enable_partition_stats();
                for (auto& state : _layers)
                    state.enable_partition_stats();
                _is_partition_stats_enabled = true;
            }
        }

        void disable_partition_stats()
        {
            BaseState::disable_partition_stats();
            for (auto& state : _layers)
                state.disable_partition_stats();
            _is_partition_stats_enabled = false;
        }

        void init_mcmc(double c, double dl)
        {
            BaseState::init_mcmc(c, dl);
            for (auto& state : _layers)
                state.init_mcmc(numeric_limits<double>::infinity(), dl);
        }
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855

        LayerState& get_layer(size_t l)
        {
            assert(l < _layers.size());
            return _layers[l];
        }

        size_t get_block(size_t l, size_t v)
        {
            return _layers[l]._b[v];
        }

        void set_block(size_t l, size_t v, size_t r)
        {
            _layers[l]._b[v] = r;
        }

        size_t get_vweight(size_t l, size_t v)
        {
            return _layers[l]._vweight[v];
        }

        void couple_state(LayeredBlockStateVirtualBase& s,
                          entropy_args_t ea)
        {
            _lcoupled_state = &s;

            entropy_args_t lea(ea);
856
            //lea.edges_dl = false;
857
858
859
860
            lea.partition_dl = false;
            for (size_t l = 0; l < _layers.size(); ++l)
                _layers[l].couple_state(s.get_layer(l), lea);

861
            lea.partition_dl = ea.partition_dl;
862
863
            lea.adjacency = false;
            lea.recs = false;
864
            lea.recs_dl = false;
865
            lea.degree_dl = false;
866
            lea.edges_dl = false;
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888

            BaseState::couple_state(s, lea);

            // assert(check_layers());
        }

        void decouple_state()
        {
            BaseState::decouple_state();
            _lcoupled_state = nullptr;
            for (auto& state : _layers)
                state.decouple_state();
        }

        void couple_state(BlockStateVirtualBase& s,
                          entropy_args_t ea)
        {
            BaseState::couple_state(s, ea);
        }

        void add_partition_node(size_t v, size_t r)
        {
889
            if (_wr[r] == 0 && BaseState::_vweight[v] > 0)
890
891
892
893
894
895
896
                _actual_B++;
            BaseState::add_partition_node(v, r);
        }

        void remove_partition_node(size_t v, size_t r)
        {
            BaseState::remove_partition_node(v, r);
897
            if (_wr[r] == 0 && BaseState::_vweight[v] > 0)
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
                _actual_B--;
        }

        size_t get_layer_node(size_t l, size_t v)
        {
            auto& ls = _vc[v];
            auto& vs = _vmap[v];

            auto pos = std::lower_bound(ls.begin(), ls.end(), l);

            if (pos == ls.end() || size_t(*pos) != l)
                return null_group;

            return *(vs.begin() + (pos - ls.begin()));
        }

        void add_layer_node(size_t l, size_t v, size_t u)
        {
            auto& ls = _vc_c[v];
            auto& vs = _vmap_c[v];

            auto pos = std::lower_bound(ls.begin(), ls.end(), l);
            assert(pos == ls.end() || size_t(*pos) != l);

            vs.insert(vs.begin() + (pos - ls.begin()), u);
            ls.insert(pos, l);

            auto& state = _layers[l];
            state.set_vertex_weight(u, 0);
        }

        void remove_layer_node(size_t l, size_t v, size_t)
        {
            auto& ls = _vc[v];
            auto& vs = _vmap[v];

            auto pos = std::lower_bound(ls.begin(), ls.end(), l);

            assert(pos != ls.end());
            assert(size_t(*pos) == l);
            //assert(u == size_t(*(vs.begin() + (pos - ls.begin()))));

            vs.erase(vs.begin() + (pos - ls.begin()));
            ls.erase(pos);
        }

        size_t get_block_map(size_t l, size_t r, bool put_new=true)
        {
            return _layers[l].get_block_map(r, put_new);
        }

        void set_vertex_weight(size_t v, int w)
        {
            if (w == 0 && BaseState::_vweight[v] > 0)
                _N--;
            if (w == 1 && BaseState::_vweight[v] == 0)
                _N++;
            BaseState::set_vertex_weight(v, w);
        }

        size_t add_block()
        {
            auto r = BaseState::add_block();
            for (size_t l = 0; l < _layers.size(); ++l)
            {
                auto& state = _layers[l];
                size_t r_u = state.add_block();
                if (_lcoupled_state != nullptr)
                    _lcoupled_state->get_layer(l).coupled_resize_vertex(r_u);
            }
            return r;
        }

        void coupled_resize_vertex(size_t v)
        {
            BaseState::coupled_resize_vertex(v);
            auto& ls = _vc_c[v];
            auto& vs = _vmap_c[v];
            for (size_t j = 0; j < ls.size(); ++j)
            {
                int l = ls[j];
                size_t u = vs[j];
                auto& state = _layers[l];
                state.coupled_resize_vertex(u);
            }
        }

        void add_edge(const GraphInterface::edge_t& e)
        {
            BaseState::add_edge(e);
        }

        void remove_edge(const GraphInterface::edge_t& e)
        {
            BaseState::remove_edge(e);
        }

995
        void add_edge_rec(const GraphInterface::edge_t& e)
996
        {
997
998
999
1000
            BaseState::add_edge_rec(e);
        }

        void remove_edge_rec(const GraphInterface::edge_t& e)
For faster browsing, not all history is shown. View entire blame