graph_blockmodel.hh 84.6 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
21
22
23
24
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_BLOCKMODEL_HH
#define GRAPH_BLOCKMODEL_HH

#include "config.h"

#include <vector>

25
#include "../support/graph_state.hh"
26
27
#include "graph_blockmodel_util.hh"

28
29
#include "openmp_lock.hh"

30
31
32
33
34
35
36
37
38
39
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;

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
template <class PMap>
auto uncheck(boost::any& amap, PMap*)
{
    return any_cast<typename PMap::checked_t&>(amap).get_unchecked();
}

template <class T, class V>
auto& uncheck(boost::any& amap, UnityPropertyMap<T,V>*)
{
    return any_cast<UnityPropertyMap<T,V>&>(amap);
}

inline simple_degs_t uncheck(boost::any& amap, simple_degs_t*)
{
    return any_cast<simple_degs_t>(amap);
}

57
58
59
60
typedef mpl::vector2<std::true_type, std::false_type> bool_tr;
typedef mpl::vector2<vcmap_t, vmap_t> vweight_tr;
typedef mpl::vector2<ecmap_t, emap_t> eweight_tr;

61
62
63
64
65
66
67
68
69
70
#ifdef GRAPH_BLOCKMODEL_RMAP_ENABLE
#    ifdef GRAPH_BLOCKMODEL_RMAP_ALL_ENABLE
typedef mpl::vector2<std::true_type, std::false_type> rmap_tr;
#    else
typedef mpl::vector1<std::true_type> rmap_tr;
#    endif
#else
typedef mpl::vector1<std::false_type> rmap_tr;
#endif

71
72
#define BLOCK_STATE_params                                                     \
    ((g, &, all_graph_views, 1))                                               \
73
74
    ((is_weighted,, bool_tr, 1))                                               \
    ((use_hash,, bool_tr, 1))                                                  \
75
    ((use_rmap,, rmap_tr, 1))                                                  \
76
    ((_abg, &, boost::any&, 0))                                                \
77
78
    ((_aeweight, &, boost::any&, 0))                                           \
    ((_avweight, &, boost::any&, 0))                                           \
79
    ((_adegs, &, boost::any&, 0))                                              \
80
81
82
83
84
    ((mrs,, emap_t, 0))                                                        \
    ((mrp,, vmap_t, 0))                                                        \
    ((mrm,, vmap_t, 0))                                                        \
    ((wr,, vmap_t, 0))                                                         \
    ((b,, vmap_t, 0))                                                          \
85
86
87
88
    ((empty_blocks, & ,std::vector<size_t>&, 0))                               \
    ((empty_pos,, vmap_t, 0))                                                  \
    ((candidate_blocks, &, std::vector<size_t>&, 0))                           \
    ((candidate_pos,, vmap_t, 0))                                              \
89
90
    ((bclabel,, vmap_t, 0))                                                    \
    ((pclabel,, vmap_t, 0))                                                    \
91
    ((bfield,, vprop_map_t<std::vector<double>>::type, 0))                     \
Tiago Peixoto's avatar
Tiago Peixoto committed
92
    ((Bfield, &, std::vector<double>&, 0))                                     \
93
94
    ((merge_map,, vmap_t, 0))                                                  \
    ((deg_corr,, bool, 0))                                                     \
95
    ((rec_types,, std::vector<int32_t>, 0))                                    \
96
97
98
99
    ((rec,, std::vector<eprop_map_t<double>::type>, 0))                        \
    ((drec,, std::vector<eprop_map_t<double>::type>, 0))                       \
    ((brec,, std::vector<eprop_map_t<double>::type>, 0))                       \
    ((bdrec,, std::vector<eprop_map_t<double>::type>, 0))                      \
100
    ((brecsum,, vprop_map_t<double>::type, 0))                                 \
101
    ((wparams, &, std::vector<std::vector<double>>&, 0))                       \
102
103
    ((recdx, &, std::vector<double>&, 0))                                      \
    ((Lrecdx, &, std::vector<double>&, 0))                                     \
104
    ((epsilon, &, std::vector<double>&, 0))                                    \
105
    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))                \
106
    ((bignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))               \
107
    ((allow_empty,, bool, 0))
108
109
110
111
112

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

template <class... Ts>
class BlockState
113
    : public BlockStateBase<Ts...>, public BlockStateVirtualBase
114
115
116
117
118
{
public:
    GET_PARAMS_USING(BlockStateBase<Ts...>, BLOCK_STATE_params)
    GET_PARAMS_TYPEDEF(Ts, BLOCK_STATE_params)

119
120
    typedef partition_stats<use_rmap_t::value> partition_stats_t;

121
122
123
124
125
126
    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()),
127
128
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
129
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
130
          _emat(_bg, rng),
131
          _neighbor_sampler(_g, _eweight),
132
          _m_entries(num_vertices(_bg))
133
    {
134
135
        _empty_blocks.clear();
        _candidate_blocks.clear();
136
        _candidate_blocks.push_back(null_group);
137
138
139
140
141
142
143
        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);
        }
144
145
146
147
        for (auto& p : _rec)
            _c_rec.push_back(p.get_checked());
        for (auto& p : _drec)
            _c_drec.push_back(p.get_checked());
148
149
150
151
152
153
154
155
156
157
158
        for (auto& p : _brec)
        {
            _c_brec.push_back(p.get_checked());
            double x = 0;
            for (auto me : edges_range(_bg))
                x += p[me];
            _recsum.push_back(x);
        }
        for (auto& p : _bdrec)
            _c_bdrec.push_back(p.get_checked());

159
        if (!_rec_types.empty())
160
        {
161
162
            _recx2.resize(_rec_types.size());
            _recdx.resize(_rec_types.size());
163
            for (auto me : edges_range(_bg))
164
            {
165
166
167
168
169
                if (_brec[0][me] > 0)
                {
                    _B_E++;
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                    {
170
                        if (_rec_types[i] == weight_type::REAL_NORMAL)
171
172
173
174
175
176
177
178
179
180
181
                        {
                            _recx2[i] += std::pow(_brec[i][me], 2);
                            if (_brec[0][me] > 1)
                                _recdx[i] += \
                                    (_bdrec[i][me] -
                                     std::pow(_brec[i][me], 2) / _brec[0][me]);
                        }
                    }
                }
                if (_brec[0][me] > 1)
                    _B_E_D++;
182
183
            }
        }
184
185
186
187
188
189
190
191
192

        _rt = weight_type::NONE;
        for (auto rt : _rec_types)
        {
            _rt = rt;
            if (rt == weight_type::REAL_NORMAL)
                break;
        }
        _dBdx.resize(_rec_types.size());
193
        _LdBdx.resize(_rec_types.size());
194
195
196
197
198
199
    }

    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()),
200
201
          _c_rec(other._c_rec),
          _c_drec(other._c_drec),
202
203
204
205
          _c_brec(other._c_brec),
          _c_bdrec(other._c_bdrec),
          _recsum(other._recsum),
          _recx2(other._recx2),
206
          _dBdx(other._dBdx),
207
          _LdBdx(other._LdBdx),
208
          _B_E(other._B_E),
209
210
          _B_E_D(other._B_E_D),
          _rt(other._rt),
211
212
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
213
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
214
          _emat(other._emat),
215
          _egroups_enabled(other._egroups_enabled),
216
          _neighbor_sampler(other._neighbor_sampler),
217
          _m_entries(num_vertices(_bg))
218
219
220
221
222
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

223
224
225
    // =========================================================================
    // State modification
    // =========================================================================
226

227
    template <class MEntries, class EFilt>
228
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
229
                          EFilt&& efilt)
230
231
    {
        auto mv_entries = [&](auto&&... args)
232
            {
233
234
235
                move_entries(v, r, nr, _b, _g, _eweight, num_vertices(_bg),
                             m_entries, std::forward<EFilt>(efilt),
                             is_loop_nop(),
236
                             std::forward<decltype(args)>(args)...);
237
            };
238

239
        if (_rt == weight_type::NONE)
240
        {
241
            mv_entries();
242
        }
243
        else
244
        {
245
246
247
248
            if (_rt == weight_type::REAL_NORMAL)
                mv_entries(_rec, _drec);
            else
                mv_entries(_rec);
249
250
        }
    }
251

252
253
254
    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
255
        get_move_entries(v, r, nr, m_entries, [](auto) {return false;});
256
    }
257
258


259
260
    template <bool Add, class EFilt>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt)
261
262
    {
        if (Add)
263
            get_move_entries(v, null_group, r, _m_entries,
264
                             std::forward<EFilt>(efilt));
265
        else
266
            get_move_entries(v, r, null_group, _m_entries,
267
                             std::forward<EFilt>(efilt));
Tiago Peixoto's avatar
Tiago Peixoto committed
268

269
        apply_delta<Add,!Add>(*this, _m_entries);
270
271
272
273
274
275
276
277
278
279

        if (Add)
            BlockState::add_partition_node(v, r);
        else
            BlockState::remove_partition_node(v, r);

        if (!_rec_types.empty() &&
            _rec_types[1] == weight_type::DELTA_T) // waiting times
        {
            if (_ignore_degrees[v] > 0)
280
            {
281
282
283
284
285
286
287
288
289
                auto dt = out_degreeS()(v, _g, _rec[1]);
                if (Add)
                    _brecsum[r] += dt;
                else
                    _brecsum[r] -= dt;
            }
        }
    }

290
    bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true)
291
    {
292
293
294
295
296
297
298
        if (_coupled_state != nullptr && is_last(v))
        {
            auto& bh = _coupled_state->get_b();
            if (bh[r] != bh[nr])
                return false;
        }

299
300
301
302
303
304
        if (allow_empty)
            return ((_bclabel[r] == _bclabel[nr]) || (_wr[nr] == 0));
        else
            return _bclabel[r] == _bclabel[nr];
    }

305
306
307
    template <class EFilt>
    void move_vertex(size_t v, size_t r, size_t nr, EFilt&& efilt)
    {
308
309
310
        if (r == nr)
            return;

311
        if (!allow_move(v, r, nr))
312
313
            throw ValueException("cannot move vertex across clabel barriers");

314
315
        get_move_entries(v, r, nr, _m_entries, std::forward<EFilt>(efilt));

316
        apply_delta<true, true>(*this, _m_entries);
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332

        BlockState::remove_partition_node(v, r);
        BlockState::add_partition_node(v, nr);

        if (!_rec_types.empty() &&
            _rec_types[1] == weight_type::DELTA_T) // waiting times
        {
            if (_ignore_degrees[v] > 0)
            {
                auto dt = out_degreeS()(v, _g, _rec[1]);
                _brecsum[r] -= dt;
                _brecsum[nr] += dt;
            }
        }
    }

333
334
    // move a vertex from its current block to block nr
    void move_vertex(size_t v, size_t r, size_t nr)
335
    {
336
337
        move_vertex(v, r, nr, [](auto&) {return false;});
    }
338

339
340
341
342
343
    void move_vertex(size_t v, size_t nr)
    {
        size_t r = _b[v];
        move_vertex(v, r, nr);
    }
344

345
346
347
348
349
350
351
352
    void propagate_delta(size_t u, size_t v,
                         std::vector<std::tuple<size_t, size_t,
                                                GraphInterface::edge_t, int,
                                                std::vector<double>>>& entries)
    {
        size_t r = _b[u];
        size_t s = _b[v];
        _m_entries.set_move(r, s, num_vertices(_bg));
353

354
        if (_rt == weight_type::NONE)
355
        {
356
            for (auto& rsd : entries)
357
            {
358
359
360
361
                _m_entries.template insert_delta<true>(_b[get<0>(rsd)],
                                                       _b[get<1>(rsd)],
                                                       get<3>(rsd));
                update_edge(get<2>(rsd));
362
            }
363
364
365
        }
        else
        {
366
            for (auto& rsd : entries)
367
            {
368
369
370
371
                recs_propagate_insert(*this, _b[get<0>(rsd)], _b[get<1>(rsd)],
                                      get<2>(rsd), get<3>(rsd), get<4>(rsd),
                                      _m_entries);
                update_edge(get<2>(rsd));
372
            }
373
        }
374
        apply_delta<true, true>(*this, _m_entries);
375
    }
376

377
    void add_edge(const GraphInterface::edge_t& e)
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
    {
        size_t r = _b[source(e, _g)];
        size_t s = _b[target(e, _g)];
        auto me = _emat.get_me(r, s);
        if (me == _emat.get_null_edge())
        {
            me = boost::add_edge(r, s, _bg).first;
            _emat.put_me(r, s, me);
            _c_mrs[me] = 0;
            for (size_t i = 0; i < _rec_types.size(); ++i)
            {
                _c_brec[i][me] = 0;
                _c_bdrec[i][me] = 0;
            }
            if (_coupled_state != nullptr)
                _coupled_state->add_edge(me);
        }
    }

    void remove_edge(const GraphInterface::edge_t& e)
    {
        size_t r = _b[source(e, _g)];
        size_t s = _b[target(e, _g)];
401
402
        auto me = _emat.get_me(r, s);
        if (me != _emat.get_null_edge() && _mrs[me] == 0)
403
404
405
406
407
        {
            _emat.remove_me(me, _bg);
            if (_coupled_state != nullptr)
                _coupled_state->remove_edge(me);
        }
408
409
        assert(e != _emat.get_null_edge());
        update_edge(e);
410
411
412
        boost::remove_edge(e, _g);
    }

413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
    void update_edge(const GraphInterface::edge_t& e)
    {
        update_edge(e, is_weighted_t());
    }

    void update_edge(const GraphInterface::edge_t& e, std::true_type)
    {
        if (e == _emat.get_null_edge())
            return;
        auto u = source(e, _g);
        auto v = target(e, _g);
        auto ew = _eweight[e];
        _neighbor_sampler.remove(u, v, e);
        if (ew > 0)
            _neighbor_sampler.insert(u, v, ew, e);
        if (u != v)
        {
            _neighbor_sampler.remove(v, u, e);
            if (ew > 0)
                _neighbor_sampler.insert(v, u, ew, e);
        }
        if (_egroups_enabled && !_egroups.empty())
        {
            _egroups.remove_edge(e, _b, _g);
            if (ew > 0)
                _egroups.insert_edge(e, ew, _b, _g);
        }
    }

    void update_edge(const GraphInterface::edge_t&, std::false_type)
    {
    }

446
    void add_edge_rec(const GraphInterface::edge_t& e)
447
    {
448
449
450
451
452
        if (_rec_types.empty())
            return;
        auto crec = _rec[0].get_checked();
        crec[e] = 1;
        for (size_t i = 1; i < _rec_types.size(); ++i)
453
        {
454
455
456
            auto drec = _drec[i].get_checked();
            drec[e] = 0;
        }
457
458
    }

459
    void remove_edge_rec(const GraphInterface::edge_t& e)
460
    {
461
462
463
        if (_rec_types.empty())
            return;
        _rec[0][e] = 0;
464
465
    }

466
467
    void update_edge_rec(const GraphInterface::edge_t& e,
                         const std::vector<double>& delta)
468
    {
469
470
471
472
        if (_rec_types.empty())
            return;

        for (size_t i = 0; i < _rec_types.size(); ++i)
473
474
475
476
        {
            if (_rec_types[i] != weight_type::REAL_NORMAL)
                continue;

477
478
479
480
            auto rec = _c_rec[i][e];
            auto d = (std::pow(rec, 2) -
                      std::pow(rec - delta[i], 2));
            _c_drec[i][e] += d;
481
482
483
        }
    }

484
485
    void remove_partition_node(size_t v, size_t r)
    {
486
487
        assert(size_t(_b[v]) == r);

488
489
490
491
492
493
494
495
496
497
498
499
500
        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
        {
            remove_element(_candidate_blocks, _candidate_pos, r);
            add_element(_empty_blocks, _empty_pos, r);

            if (_coupled_state != nullptr)
            {
                auto& hb = _coupled_state->get_b();
                _coupled_state->remove_partition_node(r, hb[r]);
                _coupled_state->set_vertex_weight(r, 0);
            }
        }

501
        _wr[r] -= _vweight[v];
502

503
        if (!_egroups.empty() && _egroups_enabled)
504
            _egroups.remove_vertex(v, _b, _g);
505
506

        if (is_partition_stats_enabled())
507
508
509
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
510
511
512
513
    }

    void add_partition_node(size_t v, size_t r)
    {
514
515
        _b[v] = r;

516
        _wr[r] += _vweight[v];
517

518
        if (!_egroups.empty() && _egroups_enabled)
519
            _egroups.add_vertex(v, _b, _eweight, _g);
520
521

        if (is_partition_stats_enabled())
522
523
            get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
                                              _eweight, _degs);
Tiago Peixoto's avatar
Tiago Peixoto committed
524

525
        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
Tiago Peixoto's avatar
Tiago Peixoto committed
526
        {
527
            remove_element(_empty_blocks, _empty_pos, r);
528
            add_element(_candidate_blocks, _candidate_pos, r);
529
530
531
532
533
534
535

            if (_coupled_state != nullptr)
            {
                auto& hb = _coupled_state->get_b();
                _coupled_state->set_vertex_weight(r, 1);
                _coupled_state->add_partition_node(r, hb[r]);
            }
536
537
538
        }
    }

539
540
541
    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
542
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt));
543
544
545
546
    }

    void remove_vertex(size_t v, size_t r)
    {
547
        remove_vertex(v, r, [](auto&) { return false; });
548
549
    }

550
551
    void remove_vertex(size_t v)
    {
552
553
        size_t r = _b[v];
        remove_vertex(v, r);
554
555
    }

556
557
    template <class Vlist>
    void remove_vertices(Vlist& vs)
558
559
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
560
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
561
562
563

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

565
566
567
568
569
570
571
572
573
574
575
        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)
576
            remove_vertex(v, _b[v],
577
                          [&](auto& e) { return eset.find(e) != eset.end(); });
578
579
580
581
582
583
584
585

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

586
            auto me = _emat.get_me(r, s);
587
588
589
590
591
592
593
594
595

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

            assert(_mrs[me] >= 0);

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

596
            for (size_t i = 0; i < _rec_types.size(); ++i)
597
            {
598
599
600
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
601
                    _bdrec[i][me] -= _drec[i][e];
602
                    [[gnu::fallthrough]];
603
604
                default:
                    _brec[i][me] -= _rec[i][e];
605
                }
606
            }
607

608
609
610
611
612
            if (_mrs[me] == 0)
            {
                _emat.remove_me(me, _bg);
                if (_coupled_state != nullptr)
                    _coupled_state->remove_edge(me);
613
614
                else
                    boost::remove_edge(me, this->_bg);
615
            }
616
617
618
        }
    }

619
620
621
622
623
624
    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

625
626
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
627
    {
628
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
629
630
    }

631
632
    void add_vertex(size_t v, size_t r)
    {
633
        add_vertex(v, r, [](auto&){ return false; });
634
635
    }

636
637
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
638
    {
639
640
641
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
        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)
663
            add_vertex(vr.first, vr.second,
664
                       [&](auto& e){ return eset.find(e) != eset.end(); });
665
666
667
668
669
670
671
672
673
674

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

675
            if (me == _emat.get_null_edge())
676
            {
677
                me = boost::add_edge(r, s, _bg).first;
678
679
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
680
681
682
683
684
                for (size_t i = 0; i < _rec_types.size(); ++i)
                {
                    _c_brec[i][me] = 0;
                    _c_bdrec[i][me] = 0;
                }
685
686
687

                if (_coupled_state != nullptr)
                    _coupled_state->add_edge(me);
688
689
690
691
692
693
694
695
696
            }

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

            auto ew = _eweight[e];

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

698
            for (size_t i = 0; i < _rec_types.size(); ++i)
699
            {
700
701
702
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
703
                    _bdrec[i][me] += _drec[i][e];
704
                    [[gnu::fallthrough]];
705
706
                default:
                    _brec[i][me] += _rec[i][e];
707
                }
708
709
710
711
            }
        }
    }

712
    void add_vertices(python::object ovs, python::object ors)
713
    {
714
715
716
717
718
        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);
    }

719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
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
791
792
793
794
795
796
797
798
799
800
801
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
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
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
    template <bool Add, bool Deplete=true>
    void modify_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                     const std::vector<double>& rec)
    {
        size_t r = _b[u];
        size_t s = _b[v];

        if (is_partition_stats_enabled())
        {
            get_partition_stats(u).remove_vertex(u, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
            if (u != v)
                get_partition_stats(v).remove_vertex(v, s, _deg_corr, _g,
                                                     _vweight, _eweight,
                                                     _degs);
        }

        auto me = _emat.get_me(r, s);
        if (Add)
        {
            if (me == _emat.get_null_edge())
            {
                me = boost::add_edge(r, s, _bg).first;
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
                for (size_t i = 0; i < _rec_types.size(); ++i)
                {
                    _c_brec[i][me] = 0;
                    _c_bdrec[i][me] = 0;
                }
            }

            if (_coupled_state == nullptr)
                _mrs[me]++;
            _mrp[r]++;
            _mrm[s]++;
        }
        else
        {
            assert(me != _emat.get_null_edge());
            if (_coupled_state == nullptr)
                _mrs[me]--;
            _mrp[r]--;
            _mrm[s]--;
        }

        // constexpr auto one = (Add) ? 1 : -1;
        // for (size_t i = 0; i < _rec_types.size(); ++i)
        // {
        //     switch (_rec_types[i])
        //     {
        //     case weight_type::REAL_NORMAL: // signed weights
        //         _bdrec[i][me] += one * std::pow(rec[i], 2);
        //         throw GraphException("Lrecdx, etc...");
        //         [[gnu::fallthrough]];
        //     default:
        //         _brec[i][me] += one * rec[i];
        //     }
        // }

        modify_edge<Add, Deplete>(u, v, e, _is_weighted);

        if (is_partition_stats_enabled())
        {
            get_partition_stats(u).add_vertex(u, r, _deg_corr, _g,
                                              _vweight, _eweight,
                                              _degs);
            if (u != v)
                get_partition_stats(v).add_vertex(v, s, _deg_corr, _g,
                                                  _vweight, _eweight,
                                                  _degs);
            get_partition_stats(u).change_E(Add ? 1 : -1); // FIXME: wrong for multiple partition stats
        }

        if (_coupled_state != nullptr)
        {
            if (Add)
                _coupled_state->add_edge(r, s, me, rec);
            else
                _coupled_state->remove_edge(r, s, me, rec);
        }
    }

    template <bool Add, bool Deplete>
    void modify_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                     std::false_type)
    {
        if (Add)
        {
            e = boost::add_edge(u, v, _g).first;
        }
        else
        {
            if (Deplete)
            {
                boost::remove_edge(e, _g);
                e = GraphInterface::edge_t();
            }
        }
    }

    template <bool Add, bool Deplete>
    void modify_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                     std::true_type)
    {
        if (Add)
        {
            if (e == GraphInterface::edge_t())
            {
                e = boost::add_edge(u, v, _g).first;
                auto c_eweight = _eweight.get_checked();
                c_eweight[e] = 1;
            }
            else
            {
                _eweight[e]++;
            }
            if (_deg_corr)
            {
                get<1>(_degs[u].front())++;
                if (graph_tool::is_directed(_g))
                    get<0>(_degs[v].front())++;
                else
                    get<1>(_degs[v].front())++;
            }
        }
        else
        {
            _eweight[e]--;
            if (_eweight[e] == 0 && Deplete)
            {
                boost::remove_edge(e, _g);
                e = GraphInterface::edge_t();
            }
            if (_deg_corr)
            {
                get<1>(_degs[u].front())--;
                if (graph_tool::is_directed(_g))
                    get<0>(_degs[v].front())--;
                else
                    get<1>(_degs[v].front())--;
            }
        }
    }

    void add_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                  const std::vector<double>& rec)
    {
        modify_edge<true, false>(u, v, e, rec);
    }

    void remove_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                     const std::vector<double>& rec)
    {
        modify_edge<false, false>(u, v, e, rec);
    }

877
878
879
880
881
882
883
884
885
886
887
888
889
890
    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;
891
892
    }

893
894
895
896
897
898
899
900
901
902
903
904
905
    void init_vertex_weight(size_t v)
    {
        init_vertex_weight(v, _vweight);
    }

    void init_vertex_weight(size_t, vcmap_t&)
    {
    }

    template <class VMap>
    void init_vertex_weight(size_t v, VMap&& vweight)
    {
        vweight.resize(num_vertices(_g));
906
        vweight[v] = 0;
907
908
    }

909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
    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);
    }

925
926
927
928
929
930
931
932
933
934
935
936
937
    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());
    }

938
939
940
941
942
943
944
945
946
947
948
949
950
951
    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>
952
    void merge_vertices(size_t u, size_t v, Emap&& ec)
953
    {
954
        merge_vertices(u, v, ec, _is_weighted);
955
956
    }

957
    template <class Emap>
958
    void merge_vertices(size_t, size_t, Emap&&, std::false_type)
959
960
961
962
963
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
964
    void merge_vertices(size_t u, size_t v, Emap&& ec, std::true_type)
965
    {
966
967
968
        if (u == v)
            return;

969
        auto eweight_c = _eweight.get_checked();
970
971
972
973
974
975
        std::vector<typename rec_t::value_type::checked_t> c_rec;
        std::vector<typename brec_t::value_type::checked_t> c_drec;
        for (auto& p : _rec)
            c_rec.push_back(p.get_checked());
        for (auto& p : _drec)
            c_drec.push_back(p.get_checked());
976
977
978
979
980
981

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

986
        size_t nrec = this->_rec_types.size();
987
        std::vector<double> ecc(nrec), decc(nrec);
988

989
990
991
992
993
994
995
        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;
996

997
998
            std::fill(ecc.begin(), ecc.end(), 0);
            std::fill(decc.begin(), decc.end(), 0);
999
            for (auto& e : es)
1000
            {
For faster browsing, not all history is shown. View entire blame