graph_blockmodel.hh 55.5 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
31
32
33
34
35
36
37
// 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_HH
#define GRAPH_BLOCKMODEL_HH

#include "config.h"

#include <vector>

#include "graph_state.hh"
#include "graph_blockmodel_util.hh"

namespace graph_tool
{
using namespace boost;
using namespace std;

typedef vprop_map_t<int32_t>::type vmap_t;
typedef eprop_map_t<int32_t>::type emap_t;
typedef UnityPropertyMap<int,GraphInterface::vertex_t> vcmap_t;
typedef UnityPropertyMap<int,GraphInterface::edge_t> ecmap_t;

38
39
40
41
42
43
template <class CMap>
CMap& uncheck(boost::any& amap, CMap*) { return any_cast<CMap&>(amap); }
vmap_t::unchecked_t uncheck(boost::any& amap, vmap_t::unchecked_t*);
emap_t::unchecked_t uncheck(boost::any& amap, emap_t::unchecked_t*);

typedef mpl::vector2<std::true_type, std::false_type> bool_tr;
44
typedef mpl::vector2<simple_degs_t, degs_map_t> degs_tr;
45
46
47
48
49
50
51
52
typedef mpl::vector2<vcmap_t, vmap_t> vweight_tr;
typedef mpl::vector2<ecmap_t, emap_t> eweight_tr;

enum weight_type
{
    NONE,
    POSITIVE,
    SIGNED,
53
54
    DISCRETE_GEOMETRIC,
    DISCRETE_POISSON,
55
56
57
    DELTA_T
};

58
59
60
61

#define BLOCK_STATE_params                                                     \
    ((g, &, all_graph_views, 1))                                               \
    ((degs,, degs_tr, 1))                                                      \
62
63
    ((is_weighted,, bool_tr, 1))                                               \
    ((use_hash,, bool_tr, 1))                                                  \
64
    ((_abg, &, boost::any&, 0))                                                \
65
66
    ((_aeweight, &, boost::any&, 0))                                           \
    ((_avweight, &, boost::any&, 0))                                           \
67
68
69
70
71
    ((mrs,, emap_t, 0))                                                        \
    ((mrp,, vmap_t, 0))                                                        \
    ((mrm,, vmap_t, 0))                                                        \
    ((wr,, vmap_t, 0))                                                         \
    ((b,, vmap_t, 0))                                                          \
72
73
74
75
    ((empty_blocks, & ,std::vector<size_t>&, 0))                               \
    ((empty_pos,, vmap_t, 0))                                                  \
    ((candidate_blocks, &, std::vector<size_t>&, 0))                           \
    ((candidate_pos,, vmap_t, 0))                                              \
76
77
78
79
    ((bclabel,, vmap_t, 0))                                                    \
    ((pclabel,, vmap_t, 0))                                                    \
    ((merge_map,, vmap_t, 0))                                                  \
    ((deg_corr,, bool, 0))                                                     \
80
81
82
83
84
85
86
87
88
89
90
91
    ((rec_type,, int, 0))                                                      \
    ((rec,, eprop_map_t<double>::type, 0))                                     \
    ((drec,, eprop_map_t<double>::type, 0))                                    \
    ((brec,, eprop_map_t<double>::type, 0))                                    \
    ((bdrec,, eprop_map_t<double>::type, 0))                                   \
    ((brecsum,, vprop_map_t<double>::type, 0))                                 \
    ((m0,, double, 0))                                                         \
    ((k0,, double, 0))                                                         \
    ((v0,, double, 0))                                                         \
    ((nu0,, double, 0))                                                        \
    ((alpha,, double, 0))                                                      \
    ((beta,, double, 0))                                                       \
92
    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))                \
93
    ((bignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))               \
94
    ((allow_empty,, bool, 0))
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

template <class... Ts>
class BlockState
    : public BlockStateBase<Ts...>
{
public:
    GET_PARAMS_USING(BlockStateBase<Ts...>, BLOCK_STATE_params)
    GET_PARAMS_TYPEDEF(Ts, BLOCK_STATE_params)

    template <class RNG, class... ATs,
              typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
    BlockState(RNG& rng, ATs&&... args)
        : BlockStateBase<Ts...>(std::forward<ATs>(args)...),
          _bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
          _c_mrs(_mrs.get_checked()),
112
113
114
115
          _c_brec(_brec.get_checked()),
          _c_bdrec(_bdrec.get_checked()),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
116
          _emat(_bg, rng),
117
          _neighbour_sampler(_g, _eweight),
118
          _m_entries(num_vertices(_bg)),
119
          _coupled_state(nullptr)
120
    {
121
122
        _empty_blocks.clear();
        _candidate_blocks.clear();
123
        _candidate_blocks.push_back(null_group);
124
125
126
127
128
129
130
        for (auto r : vertices_range(_bg))
        {
            if (_wr[r] == 0)
                add_element(_empty_blocks, _empty_pos, r);
            else
                add_element(_candidate_blocks, _candidate_pos, r);
        }
131
132
133
134
135
136
    }

    BlockState(const BlockState& other)
        : BlockStateBase<Ts...>(static_cast<const BlockStateBase<Ts...>&>(other)),
          _bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
          _c_mrs(_mrs.get_checked()),
137
138
139
140
          _c_brec(_brec.get_checked()),
          _c_bdrec(_bdrec.get_checked()),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
141
142
          _emat(other._emat),
          _neighbour_sampler(other._neighbour_sampler),
143
          _m_entries(num_vertices(_bg)),
144
          _coupled_state(nullptr)
145
146
147
148
149
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

150

151
152
153
    // =========================================================================
    // State modification
    // =========================================================================
154

155
    template <class MEntries, class EFilt>
156
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
157
                          EFilt&& efilt)
158
159
    {
        auto mv_entries = [&](auto&&... args)
160
            {
161
                move_entries(v, r, nr, _b, _g, _eweight, m_entries,
162
                             std::forward<EFilt>(efilt), is_loop_nop(),
163
                             std::forward<decltype(args)>(args)...);
164
            };
165

166
        switch (_rec_type)
167
        {
168
        case weight_type::POSITIVE: // positive weights
169
170
        case weight_type::DISCRETE_GEOMETRIC:
        case weight_type::DISCRETE_POISSON:
171
        case weight_type::DELTA_T:
172
            mv_entries(_rec);
173
174
            break;
        case weight_type::SIGNED: // positive and negative weights
175
            mv_entries(_rec, _drec);
176
177
178
179
180
            break;
        default: // no weights
            mv_entries();
        }
    }
181

182
183
184
    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
185
        get_move_entries(v, r, nr, m_entries, [](auto) {return false;});
186
    }
187
188


189
190
    template <bool Add, class EFilt>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt)
191
192
193
    {
        _m_entries.clear();
        if (Add)
194
            get_move_entries(v, null_group, r, _m_entries,
195
                             std::forward<EFilt>(efilt));
196
        else
197
            get_move_entries(v, r, null_group, _m_entries,
198
                             std::forward<EFilt>(efilt));
Tiago Peixoto's avatar
Tiago Peixoto committed
199

200
201
202
        entries_op(_m_entries, _emat,
                   [&](auto r, auto s, auto& me, auto& delta)
                   {
203
204
205
                       if (get<0>(delta) == 0) // can happen with zero-weight
                           return;             // edges

206
                       if (Add && me == _emat.get_null_edge())
207
208
209
210
211
212
213
214
215
216
217
218
                       {
                           me = add_edge(r, s, this->_bg).first;
                           _emat.put_me(r, s, me);
                           this->_c_mrs[me] = 0;
                           this->_c_brec[me] = 0;
                           this->_c_bdrec[me] = 0;
                       }

                       this->_mrs[me] += get<0>(delta);
                       this->_mrp[r] += get<0>(delta);
                       this->_mrm[s] += get<0>(delta);

219
220
221
                       assert(this->_mrs[me] >= 0);
                       assert(this->_mrp[r] >= 0);
                       assert(this->_mrm[s] >= 0);
222

223
224
225
226
227
                       switch (this->_rec_type)
                       {
                       case weight_type::SIGNED: // signed weights
                           this->_bdrec[me] += get<2>(delta);
                       case weight_type::POSITIVE: // positive weights
228
229
                       case weight_type::DISCRETE_GEOMETRIC:
                       case weight_type::DISCRETE_POISSON:
230
                       case weight_type::DELTA_T:
231
232
                           this->_brec[me] += get<1>(delta);
                       }
233
234
235
236
237

                       if (!Add && this->_mrs[me] == 0)
                       {
                           _emat.remove_me(me, this->_bg);
                       }
238
239
240
241
                   });

        if (_rec_type == weight_type::DELTA_T) // waiting times
        {
242
            if (_ignore_degrees[v] > 0)
243
            {
244
                double dt = out_degreeS()(v, _g, _rec);
245
246
247
248
                if (Add)
                    _brecsum[r] += dt;
                else
                    _brecsum[r] -= dt;
249
            }
250
251
        }

252
        if (Add)
253
        {
254
            _b[v] = r;
255
            add_partition_node(v, r);
256
        }
257
        else
258
        {
259
            remove_partition_node(v, r);
260
        }
261
262
263
264
    }

    void remove_partition_node(size_t v, size_t r)
    {
265
        _wr[r] -= _vweight[v];
266

267
        if (!_egroups.empty())
268
            _egroups.remove_vertex(v, _b, _g);
269
270

        if (is_partition_stats_enabled())
271
272
273
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
274

275
        if (_vweight[v] > 0 && _wr[r] == 0)
276
        {
277
            remove_element(_candidate_blocks, _candidate_pos, r);
278
279
280
281
282
283
            add_element(_empty_blocks, _empty_pos, r);
        }
    }

    void add_partition_node(size_t v, size_t r)
    {
284
        _wr[r] += _vweight[v];
285

286
        if (!_egroups.empty())
287
            _egroups.add_vertex(v, _b, _eweight, _g);
288
289

        if (is_partition_stats_enabled())
290
291
            get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
                                              _eweight, _degs);
Tiago Peixoto's avatar
Tiago Peixoto committed
292

293
        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
Tiago Peixoto's avatar
Tiago Peixoto committed
294
        {
295
            remove_element(_empty_blocks, _empty_pos, r);
296
            add_element(_candidate_blocks, _candidate_pos, r);
297
298
299
        }
    }

300
301
302
    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
303
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt));
304
305
306
307
308
    }

    void remove_vertex(size_t v, size_t r)
    {
        remove_vertex(v, r,  [](auto&) { return false; });
309
310
    }

311
312
    void remove_vertex(size_t v)
    {
313
314
        size_t r = _b[v];
        remove_vertex(v, r);
315
316
    }

317
318
    template <class Vlist>
    void remove_vertices(Vlist& vs)
319
320
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
321
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336

        gt_hash_set<vertex_t> vset(vs.begin(), vs.end());

        gt_hash_set<edges_t> eset;
        for (auto v : vset)
        {
            for (auto e : all_edges_range(v, _g))
            {
                auto u = (source(e, _g) == v) ? target(e, _g) : source(e, _g);
                if (vset.find(u) != vset.end())
                    eset.insert(e);
            }
        }

        for (auto v : vset)
337
            remove_vertex(v, _b[v],
338
                          [&](auto& e) { return eset.find(e) != eset.end(); });
339
340
341
342
343
344
345
346

        for (auto& e : eset)
        {
            vertex_t v = source(e, _g);
            vertex_t u = target(e, _g);
            vertex_t r = _b[v];
            vertex_t s = _b[u];

347
            auto me = _emat.get_me(r, s);
348
349
350
351
352
353
354
355
356

            auto ew = _eweight[e];
            _mrs[me] -= ew;

            assert(_mrs[me] >= 0);

            _mrp[r] -= ew;
            _mrm[s] -= ew;

357
358
359
360
361
            switch (_rec_type)
            {
            case weight_type::SIGNED: // signed weights
                _bdrec[me] -= _drec[e];
            case weight_type::POSITIVE: // positive weights
362
363
            case weight_type::DISCRETE_GEOMETRIC:
            case weight_type::DISCRETE_POISSON:
364
365
                _brec[me] -= _rec[e];
            }
366
367
368

            if (_mrs[me] == 0)
                _emat.remove_me(me, _bg);
369
370
371
        }
    }

372
373
374
375
376
377
    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

378
379
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
380
    {
381
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
382
383
    }

384
385
    void add_vertex(size_t v, size_t r)
    {
386
        add_vertex(v, r, [](auto&){ return false; });
387
388
    }

389
390
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
391
    {
392
393
394
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;

        gt_hash_map<vertex_t, size_t> vset;
        for (size_t i = 0; i < vs.size(); ++i)
            vset[vs[i]] = rs[i];

        typedef typename graph_traits<g_t>::edge_descriptor edges_t;

        gt_hash_set<edges_t> eset;
        for (auto vr : vset)
        {
            auto v = vr.first;
            for (auto e : all_edges_range(v, _g))
            {
                auto u = (source(e, _g) == v) ? target(e, _g) : source(e, _g);
                if (vset.find(u) != vset.end())
                    eset.insert(e);
            }
        }

        for (auto vr : vset)
416
            add_vertex(vr.first, vr.second,
417
                       [&](auto& e){ return eset.find(e) != eset.end(); });
418
419
420
421
422
423
424
425
426
427

        for (auto e : eset)
        {
            vertex_t v = source(e, _g);
            vertex_t u = target(e, _g);
            vertex_t r = vset[v];
            vertex_t s = vset[u];

            auto me = _emat.get_me(r, s);

428
            if (me == _emat.get_null_edge())
429
430
431
432
            {
                me = add_edge(r, s, _bg).first;
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
433
434
                _c_brec[me] = 0;
                _c_bdrec[me] = 0;
435
436
437
438
439
440
441
442
443
            }

            assert(me == _emat.get_me(r, s));

            auto ew = _eweight[e];

            _mrs[me] += ew;
            _mrp[r] += ew;
            _mrm[s] += ew;
444

445
            switch (_rec_type)
446
            {
447
448
449
            case weight_type::SIGNED: // signed weights
                _bdrec[me] += _drec[e];
            case weight_type::POSITIVE: // positive weights
450
451
            case weight_type::DISCRETE_GEOMETRIC:
            case weight_type::DISCRETE_POISSON:
452
                _brec[me] += _rec[e];
453
454
455
456
            }
        }
    }

457
    void add_vertices(python::object ovs, python::object ors)
458
    {
459
460
461
462
463
464
465
466
467
468
        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);
    }

    bool allow_move(size_t r, size_t nr, bool allow_empty = true)
    {
        if (allow_empty)
            return ((_bclabel[r] == _bclabel[nr]) || (_wr[nr] == 0));
        else
469
            return _bclabel[r] == _bclabel[nr];
470
471
    }

472
    // move a vertex from its current block to block nr
473
    void move_vertex(size_t v, size_t r, size_t nr)
474
475
476
    {
        if (r == nr)
            return;
477
478

        if (!allow_move(r, nr))
479
            throw ValueException("cannot move vertex across clabel barriers");
480

481
482
        remove_vertex(v, r, [](auto&) {return false;});
        add_vertex(v, nr, [](auto&) {return false;});
483
484
485
486
487
488
489
490
491
492
493
494
495

        if (_coupled_state != nullptr && _vweight[v] > 0)
        {
            if (_wr[r] == 0)
            {
                _coupled_state->remove_partition_node(r, _bclabel[r]);
                _coupled_state->set_vertex_weight(r, 0);
            }

            if (_wr[nr] == _vweight[v])
            {
                _coupled_state->set_vertex_weight(nr, 1);
                _coupled_state->add_partition_node(nr, _bclabel[r]);
496
                _coupled_state->_b[nr] = _bclabel[r];
497
498
499
                _bclabel[nr] = _bclabel[r];
            }
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
500
501
    }

502
503
    void move_vertex(size_t v, size_t nr)
    {
504
505
        size_t r = _b[v];
        move_vertex(v, r, nr);
506
507
    }

508
509
510
511
512
513
514
515
516
517
518
519
520
521
    void set_vertex_weight(size_t v, int w)
    {
        set_vertex_weight(v, w, _vweight);
    }

    void set_vertex_weight(size_t, int, vcmap_t&)
    {
        throw ValueException("Cannot set the weight of an unweighted state");
    }

    template <class VMap>
    void set_vertex_weight(size_t v, int w, VMap&& vweight)
    {
        vweight[v] = w;
522
523
    }

524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
    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);
    }

540
541
542
543
544
545
546
547
548
549
550
551
552
    template <class VMap>
    void set_partition(VMap&& b)
    {
        for (auto v : vertices_range(_g))
            move_vertex(v, b[v]);
    }

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

553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
    size_t virtual_remove_size(size_t v)
    {
        return _wr[_b[v]] - _vweight[v];
    }

    // merge vertex u into v
    void merge_vertices(size_t u, size_t v)
    {
        typedef typename graph_traits<g_t>::edge_descriptor edge_t;
        UnityPropertyMap<int, edge_t> dummy;
        merge_vertices(u, v, dummy);
    }

    template <class Emap>
    void merge_vertices(size_t u, size_t v, Emap& ec)
    {
        if (u == v)
            return;
571
        merge_vertices(u, v, ec, _is_weighted);
572
573
    }

574
575
    template <class Emap>
    void merge_vertices(size_t, size_t, Emap&, std::false_type)
576
577
578
579
580
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
581
    void merge_vertices(size_t u, size_t v, Emap& ec, std::true_type)
582
583
    {
        auto eweight_c = _eweight.get_checked();
584
585
        auto rec_c = _rec.get_checked();
        auto drec_c = _drec.get_checked();
586
587
588
589
590
591

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

        gt_hash_map<std::tuple<vertex_t, int>, vector<edge_t>> ns_u, ns_v;
        for(auto e : out_edges_range(u, _g))
592
            ns_u[{target(e, _g), ec[e]}].push_back(e);
593
        for(auto e : out_edges_range(v, _g))
594
            ns_v[{target(e, _g), ec[e]}].push_back(e);
595
596
597
598
599
600
601
602

        for(auto& kv : ns_u)
        {
            vertex_t t = get<0>(kv.first);
            int l = get<1>(kv.first);
            auto& es = kv.second;

            size_t w = 0;
603
            double ecc = 0, decc = 0;
604
            for (auto& e : es)
605
            {
606
                w += _eweight[e];
607
608
609
                ecc += _rec[e];
                decc += _drec[e];
            }
610
611
612
613
614
615
616
617

            if (t == u)
            {
                t = v;
                if (!is_directed::apply<g_t>::type::value)
                {
                    assert(w % 2 == 0);
                    w /= 2;
618
619
                    ecc /= 2;
                    decc /= 2;
620
621
622
                }
            }

623
            auto iter = ns_v.find({t, l});
624
625
626
627
            if (iter != ns_v.end())
            {
                auto& e = iter->second.front();
                _eweight[e] += w;
628
629
                _rec[e] += ecc;
                _drec[e] += decc;
630
631
632
633
            }
            else
            {
                auto e = add_edge(v, t, _g).first;
634
                ns_v[{t, l}].push_back(e);
635
                eweight_c[e] = w;
636
637
                rec_c[e] = ecc;
                drec_c[e] = decc;
638
639
640
641
642
643
644
645
646
647
                set_prop(ec, e, l);
            }
        }

        if (is_directed::apply<g_t>::type::value)
        {
            ns_u.clear();
            ns_v.clear();

            for(auto e : in_edges_range(v, _g))
648
                ns_v[{source(e, _g), ec[e]}].push_back(e);
649
            for(auto e : in_edges_range(u, _g))
650
                ns_u[{source(e, _g), ec[e]}].push_back(e);
651
652
653
654
655
656
657
658
659
660
661

            for(auto& kv : ns_u)
            {
                vertex_t s = get<0>(kv.first);
                int l = get<1>(kv.first);
                auto& es = kv.second;

                if (s == u)
                    continue;

                size_t w = 0;
662
                double ecc = 0, decc = 0;
663
                for (auto& e : es)
664
                {
665
                    w += _eweight[e];
666
667
668
                    ecc += _rec[e];
                    decc += _drec[e];
                }
669

670
                auto iter = ns_v.find({s, l});
671
672
673
674
                if (iter != ns_v.end())
                {
                    auto& e = iter->second.front();
                    _eweight[e] += w;
675
676
                    _rec[e] += ecc;
                    _drec[e] += decc;
677
678
679
680
                }
                else
                {
                    auto e = add_edge(s, v, _g).first;
681
                    ns_v[{s, l}].push_back(e);
682
                    eweight_c[e] = w;
683
684
                    rec_c[e] = ecc;
                    drec_c[e] = decc;
685
686
687
688
689
690
691
692
                    set_prop(ec, e, l);
                }
            }
        }

        _vweight[v] +=_vweight[u];
        _vweight[u] = 0;
        for (auto e : all_edges_range(u, _g))
693
        {
694
            _eweight[e] = 0;
695
696
697
            _rec[e] = 0;
            _drec[e] = 0;
        }
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
        clear_vertex(u, _g);
        _merge_map[u] = v;
        merge_degs(u, v, _degs);
    }

    template <class EMap, class Edge, class Val>
    void set_prop(EMap& ec, Edge& e, Val& val)
    {
        ec[e] = val;
    }

    template <class Edge, class Val>
    void set_prop(UnityPropertyMap<Val, Edge>&, Edge&, Val&)
    {
    }

    void merge_degs(size_t, size_t, const simple_degs_t&) {}

    void merge_degs(size_t u, size_t v, typename degs_map_t::unchecked_t& degs)
    {
        gt_hash_map<std::tuple<size_t, size_t>, size_t> hist;
        for (auto& kn : degs[u])
720
            hist[{get<0>(kn), get<1>(kn)}] += get<2>(kn);
721
        for (auto& kn : degs[v])
722
            hist[{get<0>(kn), get<1>(kn)}] += get<2>(kn);
723
724
725
726
727
728
729
        degs[u].clear();
        degs[v].clear();
        auto& d = degs[v];
        for (auto& kn : hist)
            d.emplace_back(get<0>(kn.first), get<1>(kn.first), kn.second);
    }

730
731
732
    // =========================================================================
    // Virtual state modification
    // =========================================================================
733
734
735
736
737
738
739

    // compute the entropy difference of a virtual move of vertex from block r
    // to nr
    template <bool exact, class MEntries>
    double virtual_move_sparse(size_t v, size_t r, size_t nr,
                               MEntries& m_entries)
    {
740
741
742
        if (r == nr)
            return 0.;

743
        double dS = entries_dS<exact>(m_entries, _mrs, _emat, _bg);
744

745
        size_t kout = out_degreeS()(v, _g, _eweight);
746
747
        size_t kin = kout;
        if (is_directed::apply<g_t>::type::value)
748
            kin = in_degreeS()(v, _g, _eweight);
749

750
        int dwr = _vweight[v];
751
752
        int dwnr = dwr;

753
754
755
756
757
        if (r == null_group && dwnr == 0)
            dwnr = 1;

        auto vt = [&](auto mrp, auto mrm, auto nr)
            {
758
                assert(mrp >= 0 && mrm >=0 && nr >= 0);
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
                if (exact)
                    return vterm_exact(mrp, mrm, nr, _deg_corr, _bg);
                else
                    return vterm(mrp, mrm, nr, _deg_corr, _bg);
            };

        if (r != null_group)
        {
            dS += vt(_mrp[r]  - kout, _mrm[r]  - kin, _wr[r]  - dwr );
            dS -= vt(_mrp[r]        , _mrm[r]       , _wr[r]        );
        }

        if (nr != null_group)
        {
            dS += vt(_mrp[nr] + kout, _mrm[nr] + kin, _wr[nr] + dwnr);
            dS -= vt(_mrp[nr]       , _mrm[nr]      , _wr[nr]       );
        }
776
777
778
779

        return dS;
    }

780
781
    template <bool exact>
    double virtual_move_sparse(size_t v, size_t r, size_t nr)
782
    {
783
        return virtual_move_sparse<exact>(v, r, nr);
784
785
    }

786
    double virtual_move_dense(size_t v, size_t r, size_t nr, bool multigraph)
787
788
789
790
791
792
793
794
795
796
    {
        if (_deg_corr)
            throw GraphException("Dense entropy for degree corrected model not implemented!");

        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;

        if (r == nr)
            return 0;

        int kin = 0, kout = 0;
797
        kout += out_degreeS()(v, _g, _eweight);
798
        if (is_directed::apply<g_t>::type::value)
799
            kin += in_degreeS()(v, _g, _eweight);
800
801
802

        vector<int> deltap(num_vertices(_bg), 0);
        int deltal = 0;
803
        for (auto e : out_edges_range(v, _g))
804
        {
805
806
            vertex_t u = target(e, _g);
            vertex_t s = _b[u];
807
            if (u == v)
808
                deltal += _eweight[e];
809
            else
810
                deltap[s] += _eweight[e];
811
812
813
814
815
        }
        if (!is_directed::apply<g_t>::type::value)
            deltal /= 2;

        vector<int> deltam(num_vertices(_bg), 0);
816
        for (auto e : in_edges_range(v, _g))
817
        {
818
            vertex_t u = source(e, _g);
819
820
            if (u == v)
                continue;
821
822
            vertex_t s = _b[u];
            deltam[s] += _eweight[e];
823
824
825
        }

        double dS = 0;
826
        int dwr = _vweight[v];
827
828
        int dwnr = dwr;

829
830
831
832
833
834
835
836
837
838
        if (r == null_group && dwnr == 0)
            dwnr = 1;

        if (nr == null_group)
        {
            std::fill(deltap.begin(), deltap.end(), 0);
            std::fill(deltam.begin(), deltam.end(), 0);
            deltal = 0;
        }

839
840
841
        double Si = 0, Sf = 0;
        for (vertex_t s = 0; s < num_vertices(_bg); ++s)
        {
842
843
            int ers = (r != null_group) ? get_beprop(r, s, _mrs, _emat) : 0;
            int enrs = (nr != null_group) ? get_beprop(nr, s, _mrs, _emat) : 0;
844
845
846
847
848

            if (!is_directed::apply<g_t>::type::value)
            {
                if (s != nr && s != r)
                {
849
850
851
852
853
854
855
856
857
858
859
                    if (r != null_group)
                    {
                        Si += eterm_dense(r,  s, ers,              _wr[r],         _wr[s], multigraph, _bg);
                        Sf += eterm_dense(r,  s, ers - deltap[s],  _wr[r] - dwr,   _wr[s], multigraph, _bg);
                    }

                    if (nr != null_group)
                    {
                        Si += eterm_dense(nr, s, enrs,             _wr[nr],        _wr[s], multigraph, _bg);
                        Sf += eterm_dense(nr, s, enrs + deltap[s], _wr[nr] + dwnr, _wr[s], multigraph, _bg);
                    }
860
861
862
863
864
865
866
867
868
869
870
871
872
                }

                if (s == r)
                {
                    Si += eterm_dense(r, r, ers,                      _wr[r],       _wr[r],       multigraph, _bg);
                    Sf += eterm_dense(r, r, ers - deltap[r] - deltal, _wr[r] - dwr, _wr[r] - dwr, multigraph, _bg);
                }

                if (s == nr)
                {
                    Si += eterm_dense(nr, nr, enrs,                       _wr[nr],        _wr[nr],        multigraph, _bg);
                    Sf += eterm_dense(nr, nr, enrs + deltap[nr] + deltal, _wr[nr] + dwnr, _wr[nr] + dwnr, multigraph, _bg);

873
874
875
876
877
                    if (r != null_group)
                    {
                        Si += eterm_dense(r, nr, ers,                          _wr[r],       _wr[nr],        multigraph, _bg);
                        Sf += eterm_dense(r, nr, ers - deltap[nr] + deltap[r], _wr[r] - dwr, _wr[nr] + dwnr, multigraph, _bg);
                    }
878
879
880
881
                }
            }
            else
            {
882
883
                int esr = (r != null_group) ? get_beprop(s, r, _mrs, _emat) : 0;
                int esnr  = (nr != null_group) ? get_beprop(s, nr, _mrs, _emat) : 0;
884
885
886

                if (s != nr && s != r)
                {
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
                    if (r != null_group)
                    {
                        Si += eterm_dense(r, s, ers            , _wr[r]      , _wr[s]      , multigraph, _bg);
                        Sf += eterm_dense(r, s, ers - deltap[s], _wr[r] - dwr, _wr[s]      , multigraph, _bg);
                        Si += eterm_dense(s, r, esr            , _wr[s]      , _wr[r]      , multigraph, _bg);
                        Sf += eterm_dense(s, r, esr - deltam[s], _wr[s]      , _wr[r] - dwr, multigraph, _bg);
                    }

                    if (nr != null_group)
                    {
                        Si += eterm_dense(nr, s, enrs            , _wr[nr]       , _wr[s]        , multigraph, _bg);
                        Sf += eterm_dense(nr, s, enrs + deltap[s], _wr[nr] + dwnr, _wr[s]        , multigraph, _bg);
                        Si += eterm_dense(s, nr, esnr            , _wr[s]        , _wr[nr]       , multigraph, _bg);
                        Sf += eterm_dense(s, nr, esnr + deltam[s], _wr[s]        , _wr[nr] + dwnr, multigraph, _bg);
                    }
902
903
904
905
906
907
908
                }

                if(s == r)
                {
                    Si += eterm_dense(r, r, ers                                  , _wr[r]      , _wr[r]      , multigraph, _bg);
                    Sf += eterm_dense(r, r, ers - deltap[r]  - deltam[r] - deltal, _wr[r] - dwr, _wr[r] - dwr, multigraph, _bg);

909
910
911
912
913
                    if (nr != null_group)
                    {
                        Si += eterm_dense(r, nr, esnr                         , _wr[r]      , _wr[nr]       , multigraph, _bg);
                        Sf += eterm_dense(r, nr, esnr - deltap[nr] + deltam[r], _wr[r] - dwr, _wr[nr] + dwnr, multigraph, _bg);
                    }
914
915
916
917
918
919
920
                }

                if(s == nr)
                {
                    Si += eterm_dense(nr, nr, esnr                                   , _wr[nr]       , _wr[nr]       , multigraph, _bg);
                    Sf += eterm_dense(nr, nr, esnr + deltap[nr] + deltam[nr] + deltal, _wr[nr] + dwnr, _wr[nr] + dwnr, multigraph, _bg);

921
922
923
924
925
                    if (r != null_group)
                    {
                        Si += eterm_dense(nr, r, esr                         , _wr[nr]       , _wr[r]      , multigraph, _bg);
                        Sf += eterm_dense(nr, r, esr + deltap[r] - deltam[nr], _wr[nr] + dwnr, _wr[r] - dwr, multigraph, _bg);
                    }
926
927
928
929
930
931
932
                }
            }
        }

        return Sf - Si + dS;
    }

933
    template <class MEntries>
934
    double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea,
935
                        MEntries& m_entries)
936
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
937
        assert(size_t(_b[v]) == r || r == null_group);
938

939
940
        if (r == nr)
            return 0;
941

942
        if (r != null_group && nr != null_group && !allow_move(r, nr))
943
944
            return std::numeric_limits<double>::infinity();

945
        m_entries.clear();
946
        get_move_entries(v, r, nr, m_entries, [](auto) { return false; });
947

948
949
950
951
952
        double dS = 0;
        if (ea.adjacency)
        {
            if (ea.dense)
            {
953
                dS = virtual_move_dense(v, r, nr, ea.multigraph);
954
955
956
957
958
959
960
961
962
            }
            else
            {
                if (ea.exact)
                    dS = virtual_move_sparse<true>(v, r, nr, m_entries);
                else
                    dS = virtual_move_sparse<false>(v, r, nr, m_entries);
            }
        }
963

964
        if (ea.partition_dl || ea.degree_dl || ea.edges_dl)
965
966
967
        {
            enable_partition_stats();
            auto& ps = get_partition_stats(v);
968
            if (ea.partition_dl)
969
                dS += ps.get_delta_partition_dl(v, r, nr, _vweight);
970
            if (_deg_corr && ea.degree_dl)
971
972
                dS += ps.get_delta_deg_dl(v, r, nr, _vweight, _eweight,
                                          _degs, _g, ea.degree_dl_kind);
973
            if (ea.edges_dl)
974
975
976
977
            {
                size_t actual_B = 0;
                for (auto& ps : _partition_stats)
                    actual_B += ps.get_actual_B();
978
979
                dS += ps.get_delta_edges_dl(v, r, nr, _vweight, actual_B,
                                            _g);
980
            }
981
982
        }

983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
        auto positive_entries_op = [&](auto&& w_log_P)
            {
                entries_op(m_entries, this->_emat,
                           [&](auto, auto, auto& me, auto& delta)
                           {
                               size_t ers = 0;
                               double xrs = 0;
                               if (me != _emat.get_null_edge())
                               {
                                   ers = this->_mrs[me];
                                   xrs = this->_brec[me];
                               }
                               auto d = get<0>(delta);
                               auto dx = get<1>(delta);
                               dS -= -w_log_P(ers, xrs);
                               dS += -w_log_P(ers + d, xrs + dx);
                           });
            };

1002
        switch (_rec_type)
1003
        {
1004
        case weight_type::POSITIVE: // positive weights
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
            positive_entries_op([&](auto N, auto x)
                                { return positive_w_log_P(N, x,
                                                          this->_alpha,
                                                          this->_beta); });
            break;
        case weight_type::DISCRETE_GEOMETRIC:
            positive_entries_op([&](auto N, auto x)
                                { return geometric_w_log_P(N, x,
                                                           this->_alpha,
                                                           this->_beta); });
            break;
        case weight_type::DISCRETE_POISSON:
            positive_entries_op([&](auto N, auto x)
                                { return poisson_w_log_P(N, x,
                                                         this->_alpha,
                                                         this->_beta); });
1021
1022
            break;
        case weight_type::SIGNED: // positive and negative weights
1023
            entries_op(m_entries, _emat,
1024
1025
1026
1027
                       [&](auto, auto, auto& me, auto& delta)
                       {
                           size_t ers = 0;
                           double xrs = 0, x2rs = 0;
1028
                           if (me != _emat.get_null_edge())
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
                           {
                               ers = this->_mrs[me];
                               xrs = this->_brec[me];
                               x2rs = this->_bdrec[me];
                           }
                           auto d = get<0>(delta);
                           auto dx = get<1>(delta);
                           auto dx2 = get<2>(delta);
                           auto sigma1 = x2rs - xrs * (xrs / ers);
                           auto sigma2 = x2rs + dx2 - (xrs + dx) * ((xrs + dx) / (ers + d));
                           dS -= -signed_w_log_P(ers, xrs, sigma1,
                                                 this->_m0,
                                                 this->_k0,
                                                 this->_v0,
                                                 this->_nu0);
                           dS += -signed_w_log_P(ers + d, xrs + dx,
                                                 sigma2,
                                                 this->_m0,
                                                 this->_k0,
                                                 this->_v0,
                                                 this->_nu0);
                       });
            break;
        case weight_type::DELTA_T: // waiting times
1053
            if ((r != nr) && _ignore_degrees[v] > 0)
Tiago Peixoto's avatar
Tiago Peixoto committed
1054
            {
1055
1056
                double dt = out_degreeS()(v, _g, _rec);
                int k = out_degreeS()(v, _g, _eweight);
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
                if (r != null_group)
                {
                    dS -= -positive_w_log_P(_mrp[r], _brecsum[r],
                                            _alpha, _beta);
                    dS += -positive_w_log_P(_mrp[r] - k, _brecsum[r] - dt,
                                            _alpha, _beta);
                }
                if (nr != null_group)
                {
                    dS -= -positive_w_log_P(_mrp[nr], _brecsum[nr],
                                            _alpha, _beta);
                    dS += -positive_w_log_P(_mrp[nr] + k, _brecsum[nr] + dt,
                                            _alpha, _beta);
                }
Tiago Peixoto's avatar
Tiago Peixoto committed
1071
            }
1072
            break;
1073
1074
        }

1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
        if (_coupled_state != nullptr && _vweight[v] > 0)
        {
            assert(r == null_group || nr == null_group || allow_move(r, nr));
            bool r_vacate = (r != null_group) && (_wr[r] == _vweight[v]);
            bool nr_occupy = (nr != null_group) && (_wr[nr] == 0);
            if (r_vacate != nr_occupy)
            {
                if (r_vacate)
                {
                    dS += _coupled_state->virtual_move(r,
                                                       _bclabel[r],
                                                       null_group,
                                                       _coupled_entropy_args);
                }

                if (nr_occupy)
                {
                    assert(_coupled_state->_vweight[nr] == 0);
                    dS += _coupled_state->virtual_move(nr,
                                                       null_group,
                                                       _bclabel[r],
                                                       _coupled_entropy_args);
                }
            }
        }
1100
1101
1102
        return dS;
    }

1103
    double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea)
1104
    {
1105
        return virtual_move(v, r, nr, ea, _m_entries);
1106
1107
    }

1108
    double get_delta_partition_dl(size_t v, size_t r, size_t nr)
1109
1110
1111
    {
        enable_partition_stats();
        auto& ps = get_partition_stats(v);
1112
        return ps.get_delta_partition_dl(v, r, nr, _vweight);
1113
1114
    }

1115
1116
1117
1118
    // =========================================================================
    // Move proposals
    // =========================================================================

1119
1120
    // Sample node placement
    template <class RNG>
1121
    size_t sample_block(size_t v, double c, RNG& rng)
1122
1123
    {
        // attempt random block
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
        size_t s;
        if (_empty_blocks.empty())
        {
            s = uniform_sample(_candidate_blocks.begin() + 1,
                               _candidate_blocks.end(),
                               rng);
        }
        else
        {
            s = uniform_sample(_candidate_blocks, rng);
            if (s == null_group)
                s = uniform_sample(_empty_blocks, rng);
        }
1137

1138
        if (!std::isinf(c) && !_neighbour_sampler.empty(v))
1139
        {
1140
            auto u = _neighbour_sampler.sample(v, rng);
1141
1142
1143
1144
            size_t t = _b[u];
            double p_rand = 0;
            if (c > 0)
            {
1145
1146
                size_t B = (_empty_blocks.empty()) ?
                    _candidate_blocks.size() - 1 : _candidate_blocks.size();
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
                if (is_directed::apply<g_t>::type::value)
                    p_rand = c * B / double(_mrp[t] + _mrm[t] + c * B);
                else
                    p_rand = c * B / double(_mrp[t] + c * B);
            }

            typedef std::uniform_real_distribution<> rdist_t;
            if (c == 0 || rdist_t()(rng) >= p_rand)
            {
                if (_egroups.empty())
                    _egroups.init(_b, _eweight, _g, _bg);
                const auto& e = _egroups.sample_edge(t, rng);
                s = _b[target(e, _g)];
                if (s == t)
                    s = _b[source(e, _g)];
1162
                else
1163
                    assert(size_t(_b[source(e, _g)]) == t);
1164
1165
1166
1167
1168
1169
            }
        }

        return s;
    }