graph_blockmodel.hh 58.6 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
120
          _coupled_state(nullptr),
          _gstate(this)
121
    {
122
123
        _empty_blocks.clear();
        _candidate_blocks.clear();
124
        _candidate_blocks.push_back(null_group);
125
126
127
128
129
130
131
        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);
        }
132
133
134
135
136
137
    }

    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()),
138
139
140
141
          _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())),
142
143
          _emat(other._emat),
          _neighbour_sampler(other._neighbour_sampler),
144
          _m_entries(num_vertices(_bg)),
145
146
          _coupled_state(nullptr),
          _gstate(other._gstate == &other ? this : other._gstate)
147
148
149
150
151
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

152

153
154
155
    // =========================================================================
    // State modification
    // =========================================================================
156

157
    template <class MEntries, class EFilt, class GetB>
158
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
159
                          EFilt&& efilt, GetB&& get_b)
160
161
162
    {
        auto& gs = *_gstate;
        auto mv_entries = [&](auto&&... args)
163
            {
164
165
166
                move_entries(v, r, nr, std::forward<GetB>(get_b), gs._g,
                             gs._eweight, m_entries,
                             std::forward<EFilt>(efilt), is_loop_nop(),
167
                             std::forward<decltype(args)>(args)...);
168
            };
169

170
        switch (_rec_type)
171
        {
172
        case weight_type::POSITIVE: // positive weights
173
174
        case weight_type::DISCRETE_GEOMETRIC:
        case weight_type::DISCRETE_POISSON:
175
        case weight_type::DELTA_T:
176
177
178
179
180
181
182
183
184
            mv_entries(gs._rec);
            break;
        case weight_type::SIGNED: // positive and negative weights
            mv_entries(gs._rec, gs._drec);
            break;
        default: // no weights
            mv_entries();
        }
    }
185

186
187
188
189
190
    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
        get_move_entries(v, r, nr, m_entries,
                         [](auto) { return false; },
191
                         [&](auto u) -> size_t { return this->_b[u]; });
192
    }
193
194


195
196
197
    template <bool Add, class EFilt, class GetB, class BEop>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt, GetB&& get_b,
                       BEop&& beop)
198
199
200
    {
        _m_entries.clear();
        if (Add)
201
202
203
            get_move_entries(v, null_group, r, _m_entries,
                             std::forward<EFilt>(efilt),
                             std::forward<GetB>(get_b));
204
        else
205
206
207
            get_move_entries(v, r, null_group, _m_entries,
                             std::forward<EFilt>(efilt),
                             std::forward<GetB>(get_b));
Tiago Peixoto's avatar
Tiago Peixoto committed
208

209
210
211
        entries_op(_m_entries, _emat,
                   [&](auto r, auto s, auto& me, auto& delta)
                   {
212
213
                       beop(false, me);

214
215
216
217
218
219
220
221
222
223
224
225
226
                       if (Add && me == this->_emat.get_null_edge())
                       {
                           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);

227
228
229
                       assert(this->_mrs[me] >= 0);
                       assert(this->_mrp[r] >= 0);
                       assert(this->_mrm[s] >= 0);
230

231
232
233
234
235
                       switch (this->_rec_type)
                       {
                       case weight_type::SIGNED: // signed weights
                           this->_bdrec[me] += get<2>(delta);
                       case weight_type::POSITIVE: // positive weights
236
237
                       case weight_type::DISCRETE_GEOMETRIC:
                       case weight_type::DISCRETE_POISSON:
238
                       case weight_type::DELTA_T:
239
240
                           this->_brec[me] += get<1>(delta);
                       }
241
242

                       beop(true, me);
243
244
245
246
                   });

        if (_rec_type == weight_type::DELTA_T) // waiting times
        {
247
248
            auto& gs = *_gstate;
            if (gs._ignore_degrees[v] > 0)
249
            {
250
                double dt = out_degreeS()(v, gs._g, gs._rec);
251
252
253
254
                if (Add)
                    _brecsum[r] += dt;
                else
                    _brecsum[r] -= dt;
255
            }
256
257
        }

258
        if (Add)
259
        {
260
            auto&& b_v = get_b(v);
261
            b_v = r;
262
            add_partition_node(v, r);
263
        }
264
        else
265
        {
266
            remove_partition_node(v, r);
267
        }
268
269
270
271
    }

    void remove_partition_node(size_t v, size_t r)
    {
272
273
274
        auto& gs = *_gstate;

        _wr[r] -= gs._vweight[v];
275

276
277
        if (!_egroups.empty() && _gstate == this)
            _egroups.remove_vertex(v, _b, _g);
278
279

        if (is_partition_stats_enabled())
280
281
282
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, gs._g,
                                                 gs._vweight, gs._eweight,
                                                 gs._degs);
283

284
        if (gs._vweight[v] > 0 && _wr[r] == 0)
285
        {
286
            remove_element(_candidate_blocks, _candidate_pos, r);
287
288
289
290
291
292
            add_element(_empty_blocks, _empty_pos, r);
        }
    }

    void add_partition_node(size_t v, size_t r)
    {
293
294
295
        auto& gs = *_gstate;

        _wr[r] += gs._vweight[v];
296

297
298
        if (!_egroups.empty() && _gstate == this)
            _egroups.add_vertex(v, _b, _eweight, _g);
299
300

        if (is_partition_stats_enabled())
301
302
            get_partition_stats(v).add_vertex(v, r, _deg_corr, gs._g, gs._vweight,
                                              gs._eweight, gs._degs);
Tiago Peixoto's avatar
Tiago Peixoto committed
303

304
        if (gs._vweight[v] > 0 && _wr[r] == gs._vweight[v])
Tiago Peixoto's avatar
Tiago Peixoto committed
305
        {
306
            remove_element(_empty_blocks, _empty_pos, r);
307
            add_element(_candidate_blocks, _candidate_pos, r);
308
309
310
        }
    }

311
312
313
    template <class EFilt, class GetB, class BEop>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt, GetB&& get_b,
                       BEop&& beop)
314
    {
315
316
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt),
                             std::forward<GetB>(get_b), beop);
317
318
319
320
321
    }

    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
322
        remove_vertex(v, r, std::forward<EFilt>(efilt),
323
324
                      [&](auto u) -> auto& { return this->_b[u]; },
                      [](bool, const auto &){});
325
326
327
328
329
    }

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

332
333
    void remove_vertex(size_t v)
    {
334
335
        size_t r = _b[v];
        remove_vertex(v, r);
336
337
    }

338
339
    template <class Vlist>
    void remove_vertices(Vlist& vs)
340
341
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
342
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357

        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)
358
            remove_vertex(v, _b[v],
359
                          [&](auto& e) { return eset.find(e) != eset.end(); });
360
361
362
363
364
365
366
367

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

368
            auto me = _emat.get_me(r, s);
369
370
371
372
373
374
375
376
377

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

            assert(_mrs[me] >= 0);

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

378
379
380
381
382
            switch (_rec_type)
            {
            case weight_type::SIGNED: // signed weights
                _bdrec[me] -= _drec[e];
            case weight_type::POSITIVE: // positive weights
383
384
            case weight_type::DISCRETE_GEOMETRIC:
            case weight_type::DISCRETE_POISSON:
385
386
                _brec[me] -= _rec[e];
            }
387
388
389
        }
    }

390
391
392
393
394
395
    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

396
397
    template <class EFilt, class GetB, class BEop>
    void add_vertex(size_t v, size_t r, EFilt&& efilt, GetB&& get_b,
398
                    BEop&& beop)
399
    {
400
401
402
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt),
                            std::forward<GetB>(get_b), std::forward<BEop>(beop));
        //assert(size_t(get_b(v)) == r);
403
404
    }

405
406
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
407
    {
408
        add_vertex(v, r, std::forward<EFilt>(efilt),
409
410
                   [&](auto u) -> auto& { return this->_b[u]; },
                   [](bool, const auto&){});
411
412
    }

413
414
    void add_vertex(size_t v, size_t r)
    {
415
        add_vertex(v, r, [](auto&){ return false; });
416
417
    }

418
419
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
420
    {
421
422
423
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
        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)
445
            add_vertex(vr.first, vr.second,
446
                       [&](auto& e){ return eset.find(e) != eset.end(); });
447
448
449
450
451
452
453
454
455
456

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

457
            if (me == _emat.get_null_edge())
458
459
460
461
            {
                me = add_edge(r, s, _bg).first;
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
462
463
                _c_brec[me] = 0;
                _c_bdrec[me] = 0;
464
465
466
467
468
469
470
471
472
            }

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

            auto ew = _eweight[e];

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

474
            switch (_rec_type)
475
            {
476
477
478
            case weight_type::SIGNED: // signed weights
                _bdrec[me] += _drec[e];
            case weight_type::POSITIVE: // positive weights
479
480
            case weight_type::DISCRETE_GEOMETRIC:
            case weight_type::DISCRETE_POISSON:
481
                _brec[me] += _rec[e];
482
483
484
485
            }
        }
    }

486
    void add_vertices(python::object ovs, python::object ors)
487
    {
488
489
490
491
492
493
494
495
496
497
        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
498
            return _bclabel[r] == _bclabel[nr];
499
500
    }

501
    // move a vertex from its current block to block nr
502
503
    template <class GetB, class BEop>
    void move_vertex(size_t v, size_t r, size_t nr, GetB&& get_b, BEop&& beop)
504
505
506
    {
        if (r == nr)
            return;
507
508

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

511
512
513
514
        remove_vertex(v, r, [](auto&) {return false;}, std::forward<GetB>(get_b),
                      std::forward<BEop>(beop));
        add_vertex(v, nr, [](auto&) {return false;}, std::forward<GetB>(get_b),
                   std::forward<BEop>(beop));
515
516
517
518
519
520
521
522
523
524
525
526
527

        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]);
528
                _coupled_state->_b[nr] = _bclabel[r];
529
530
531
                _bclabel[nr] = _bclabel[r];
            }
        }
532
        //assert(size_t(get_b(v)) == nr);
533
534
    }

535
536
    template <class GetB, class BEop>
    void move_vertex(size_t v, size_t nr, GetB&& get_b, BEop&& beop)
Tiago Peixoto's avatar
Tiago Peixoto committed
537
538
    {
        size_t r = get_b(v);
539
540
        move_vertex(v, r, nr, std::forward<GetB>(get_b),
                    std::forward<BEop>(beop));
Tiago Peixoto's avatar
Tiago Peixoto committed
541
542
    }

543
544
    void move_vertex(size_t v, size_t nr)
    {
545
546
547
        move_vertex(v, nr,
                    [&](auto u) -> auto& { return this->_b[u]; },
                    [](bool, auto&){});
548
549
    }

550
551
552
553
554
555
556
557
558
559
560
561
562
563
    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;
564
565
    }

566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
    template <class Edge>
    void set_edge_weight(Edge&& e, int w)
    {
        set_edge_weight(e, w, _eweight);
    }

    template <class Edge>
    void set_edge_weight(Edge&&, int, ecmap_t&)
    {
        throw ValueException("Cannot set the weight of an unweighted state");
    }

    template <class Edge, class EMap>
    void set_edge_weight(Edge&& e, int w, EMap&& eweight)
    {
        auto ew_c = eweight.get_checked();
        ew_c[e] = w;
    }

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

601
602
603
604
605
606
607
608
609
610
611
612
613
    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());
    }

614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
    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;
632
        merge_vertices(u, v, ec, _is_weighted);
633
634
    }

635
636
    template <class Emap>
    void merge_vertices(size_t, size_t, Emap&, std::false_type)
637
638
639
640
641
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
642
    void merge_vertices(size_t u, size_t v, Emap& ec, std::true_type)
643
644
    {
        auto eweight_c = _eweight.get_checked();
645
646
        auto rec_c = _rec.get_checked();
        auto drec_c = _drec.get_checked();
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663

        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))
            ns_u[std::make_tuple(target(e, _g), ec[e])].push_back(e);
        for(auto e : out_edges_range(v, _g))
            ns_v[std::make_tuple(target(e, _g), ec[e])].push_back(e);

        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;
664
            double ecc = 0, decc = 0;
665
            for (auto& e : es)
666
            {
667
                w += _eweight[e];
668
669
670
                ecc += _rec[e];
                decc += _drec[e];
            }
671
672
673
674
675
676
677
678

            if (t == u)
            {
                t = v;
                if (!is_directed::apply<g_t>::type::value)
                {
                    assert(w % 2 == 0);
                    w /= 2;
679
680
                    ecc /= 2;
                    decc /= 2;
681
682
683
684
685
686
687
688
                }
            }

            auto iter = ns_v.find(std::make_tuple(t, l));
            if (iter != ns_v.end())
            {
                auto& e = iter->second.front();
                _eweight[e] += w;
689
690
                _rec[e] += ecc;
                _drec[e] += decc;
691
692
693
694
695
696
            }
            else
            {
                auto e = add_edge(v, t, _g).first;
                ns_v[std::make_tuple(t, l)].push_back(e);
                eweight_c[e] = w;
697
698
                rec_c[e] = ecc;
                drec_c[e] = decc;
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
                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))
                ns_v[std::make_tuple(source(e, _g), ec[e])].push_back(e);
            for(auto e : in_edges_range(u, _g))
                ns_u[std::make_tuple(source(e, _g), ec[e])].push_back(e);

            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;
723
                double ecc = 0, decc = 0;
724
                for (auto& e : es)
725
                {
726
                    w += _eweight[e];
727
728
729
                    ecc += _rec[e];
                    decc += _drec[e];
                }
730
731
732
733
734
735

                auto iter = ns_v.find(std::make_tuple(s, l));
                if (iter != ns_v.end())
                {
                    auto& e = iter->second.front();
                    _eweight[e] += w;
736
737
                    _rec[e] += ecc;
                    _drec[e] += decc;
738
739
740
741
742
743
                }
                else
                {
                    auto e = add_edge(s, v, _g).first;
                    ns_v[std::make_tuple(s, l)].push_back(e);
                    eweight_c[e] = w;
744
745
                    rec_c[e] = ecc;
                    drec_c[e] = decc;
746
747
748
749
750
751
752
753
                    set_prop(ec, e, l);
                }
            }
        }

        _vweight[v] +=_vweight[u];
        _vweight[u] = 0;
        for (auto e : all_edges_range(u, _g))
754
        {
755
            _eweight[e] = 0;
756
757
758
            _rec[e] = 0;
            _drec[e] = 0;
        }
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
        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])
            hist[make_tuple(get<0>(kn), get<1>(kn))] += get<2>(kn);
        for (auto& kn : degs[v])
            hist[make_tuple(get<0>(kn), get<1>(kn))] += get<2>(kn);
        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);
    }

791
792
793
    // =========================================================================
    // Virtual state modification
    // =========================================================================
794
795
796
797
798
799
800

    // 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)
    {
801
802
803
        if (r == nr)
            return 0.;

804
        double dS = entries_dS<exact>(m_entries, _mrs, _emat, _bg);
805

806
807
        auto& gs = *_gstate;
        size_t kout = out_degreeS()(v, gs._g, gs._eweight);
808
809
        size_t kin = kout;
        if (is_directed::apply<g_t>::type::value)
810
            kin = in_degreeS()(v, gs._g, gs._eweight);
811

812
        int dwr = gs._vweight[v];
813
814
        int dwnr = dwr;

815
816
817
818
819
        if (r == null_group && dwnr == 0)
            dwnr = 1;

        auto vt = [&](auto mrp, auto mrm, auto nr)
            {
820
                assert(mrp >= 0 && mrm >=0 && nr >= 0);
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
                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]       );
        }
838
839
840
841

        return dS;
    }

842
843
    template <bool exact>
    double virtual_move_sparse(size_t v, size_t r, size_t nr)
844
    {
845
        return virtual_move_sparse<exact>(v, r, nr);
846
847
    }

848
    template <class GetB>
849
    double virtual_move_dense(size_t v, size_t r, size_t nr, bool multigraph,
850
                              GetB&& get_b)
851
852
853
854
855
856
857
858
859
    {
        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;

860
        auto& gs = *_gstate;
861
        int kin = 0, kout = 0;
862
        kout += out_degreeS()(v, gs._g, gs._eweight);
863
        if (is_directed::apply<g_t>::type::value)
864
            kin += in_degreeS()(v, gs._g, gs._eweight);
865
866
867

        vector<int> deltap(num_vertices(_bg), 0);
        int deltal = 0;
868
        for (auto e : out_edges_range(v, gs._g))
869
        {
870
871
            vertex_t u = target(e, gs._g);
            vertex_t s = get_b(u);
872
            if (u == v)
873
                deltal += gs._eweight[e];
874
            else
875
                deltap[s] += gs._eweight[e];
876
877
878
879
880
        }
        if (!is_directed::apply<g_t>::type::value)
            deltal /= 2;

        vector<int> deltam(num_vertices(_bg), 0);
881
        for (auto e : in_edges_range(v, gs._g))
882
        {
883
            vertex_t u = source(e, gs._g);
884
885
            if (u == v)
                continue;
886
887
            vertex_t s = get_b(u);
            deltam[s] += gs._eweight[e];
888
889
890
        }

        double dS = 0;
891
        int dwr = gs._vweight[v];
892
893
        int dwnr = dwr;

894
895
896
897
898
899
900
901
902
903
        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;
        }

904
905
906
        double Si = 0, Sf = 0;
        for (vertex_t s = 0; s < num_vertices(_bg); ++s)
        {
907
908
            int ers = (r != null_group) ? get_beprop(r, s, _mrs, _emat) : 0;
            int enrs = (nr != null_group) ? get_beprop(nr, s, _mrs, _emat) : 0;
909
910
911
912
913

            if (!is_directed::apply<g_t>::type::value)
            {
                if (s != nr && s != r)
                {
914
915
916
917
918
919
920
921
922
923
924
                    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);
                    }
925
926
927
928
929
930
931
932
933
934
935
936
937
                }

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

938
939
940
941
942
                    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);
                    }
943
944
945
946
                }
            }
            else
            {
947
948
                int esr = (r != null_group) ? get_beprop(s, r, _mrs, _emat) : 0;
                int esnr  = (nr != null_group) ? get_beprop(s, nr, _mrs, _emat) : 0;
949
950
951

                if (s != nr && s != r)
                {
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
                    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);
                    }
967
968
969
970
971
972
973
                }

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

974
975
976
977
978
                    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);
                    }
979
980
981
982
983
984
985
                }

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

986
987
988
989
990
                    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);
                    }
991
992
993
994
995
996
997
                }
            }
        }

        return Sf - Si + dS;
    }

998
    template <class MEntries, class GetB>
999
    double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea,
1000
                        MEntries& m_entries, GetB&& get_b)
For faster browsing, not all history is shown. View entire blame