graph_blockmodel.hh 84.4 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2019 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_BLOCKMODEL_HH
#define GRAPH_BLOCKMODEL_HH

#include "config.h"

#include <vector>

25
26
27
28
#ifdef __clang__
#include <boost/algorithm/minmax_element.hpp>
#endif

29
#include "../support/graph_state.hh"
30
31
#include "graph_blockmodel_util.hh"

32
33
#include "openmp_lock.hh"

34
35
36
37
38
39
40
41
42
43
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;

44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
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);
}

61
62
63
64
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;

65
66
67
68
69
70
71
72
73
74
#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

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

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

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

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

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

160
        if (!_rec_types.empty())
161
        {
162
163
            _recx2.resize(_rec_types.size());
            _recdx.resize(_rec_types.size());
164
            for (auto me : edges_range(_bg))
165
            {
166
167
168
169
170
                if (_brec[0][me] > 0)
                {
                    _B_E++;
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                    {
171
                        if (_rec_types[i] == weight_type::REAL_NORMAL)
172
173
174
175
176
177
178
179
180
181
182
                        {
                            _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++;
183
184
            }
        }
185
186
187
188
189
190
191
192
193

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

        _N = 0;
        for (auto v : vertices_range(_g))
            _N += _vweight[v];
199
200
201
202
203
204
    }

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

229
230
231
    // =========================================================================
    // State modification
    // =========================================================================
232

233
    template <class MEntries, class EFilt>
234
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
235
                          EFilt&& efilt)
236
237
    {
        auto mv_entries = [&](auto&&... args)
238
            {
239
240
241
                move_entries(v, r, nr, _b, _g, _eweight, num_vertices(_bg),
                             m_entries, std::forward<EFilt>(efilt),
                             is_loop_nop(),
242
                             std::forward<decltype(args)>(args)...);
243
            };
244

245
        if (_rt == weight_type::NONE)
246
        {
247
            mv_entries();
248
        }
249
        else
250
        {
251
252
253
254
            if (_rt == weight_type::REAL_NORMAL)
                mv_entries(_rec, _drec);
            else
                mv_entries(_rec);
255
256
        }
    }
257

258
259
260
    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
261
        get_move_entries(v, r, nr, m_entries, [](auto) {return false;});
262
    }
263
264


265
266
    template <bool Add, class EFilt>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt)
267
268
    {
        if (Add)
269
            get_move_entries(v, null_group, r, _m_entries,
270
                             std::forward<EFilt>(efilt));
271
        else
272
            get_move_entries(v, r, null_group, _m_entries,
273
                             std::forward<EFilt>(efilt));
Tiago Peixoto's avatar
Tiago Peixoto committed
274

275
        apply_delta<Add,!Add>(*this, _m_entries);
276
277
278
279
280
281
282

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

283
    bool allow_move(size_t r, size_t nr)
284
    {
285
        if (_coupled_state != nullptr)
286
        {
287
288
289
290
            auto& hb = _coupled_state->get_b();
            auto rr = hb[r];
            auto ss = hb[nr];
            if (rr != ss && !_coupled_state->allow_move(rr, ss))
291
292
                return false;
        }
293
        return _bclabel[r] == _bclabel[nr];
294
295
    }

296
297
298
    template <class EFilt>
    void move_vertex(size_t v, size_t r, size_t nr, EFilt&& efilt)
    {
299
300
301
        if (r == nr)
            return;

302
        if (!allow_move(r, nr))
303
304
            throw ValueException("cannot move vertex across clabel barriers");

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

307
        apply_delta<true, true>(*this, _m_entries);
308
309
310
311
312

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

313
314
    // move a vertex from its current block to block nr
    void move_vertex(size_t v, size_t r, size_t nr)
315
    {
316
317
        move_vertex(v, r, nr, [](auto&) {return false;});
    }
318

319
320
321
322
323
    void move_vertex(size_t v, size_t nr)
    {
        size_t r = _b[v];
        move_vertex(v, r, nr);
    }
324

325
326
327
328
329
330
331
332
    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));
333

334
        if (_rt == weight_type::NONE)
335
        {
336
            for (auto& rsd : entries)
337
            {
338
339
340
                _m_entries.template insert_delta<true>(_b[get<0>(rsd)],
                                                       _b[get<1>(rsd)],
                                                       get<3>(rsd));
341
            }
342
343
344
        }
        else
        {
345
            for (auto& rsd : entries)
346
            {
347
348
349
                recs_propagate_insert(*this, _b[get<0>(rsd)], _b[get<1>(rsd)],
                                      get<2>(rsd), get<3>(rsd), get<4>(rsd),
                                      _m_entries);
350
            }
351
        }
352
        apply_delta<true, true>(*this, _m_entries);
353
    }
354

355
    void add_edge(const GraphInterface::edge_t& e)
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    {
        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)];
379
380
        auto me = _emat.get_me(r, s);
        if (me != _emat.get_null_edge() && _mrs[me] == 0)
381
382
383
384
385
        {
            _emat.remove_me(me, _bg);
            if (_coupled_state != nullptr)
                _coupled_state->remove_edge(me);
        }
386
        assert(e != _emat.get_null_edge());
387
388
389
390
        boost::remove_edge(e, _g);
    }

    void add_edge_rec(const GraphInterface::edge_t& e)
391
    {
392
393
394
395
396
        if (_rec_types.empty())
            return;
        auto crec = _rec[0].get_checked();
        crec[e] = 1;
        for (size_t i = 1; i < _rec_types.size(); ++i)
397
        {
398
399
400
            auto drec = _drec[i].get_checked();
            drec[e] = 0;
        }
401
402
    }

403
    void remove_edge_rec(const GraphInterface::edge_t& e)
404
    {
405
406
407
        if (_rec_types.empty())
            return;
        _rec[0][e] = 0;
408
409
    }

410
411
    void update_edge_rec(const GraphInterface::edge_t& e,
                         const std::vector<double>& delta)
412
    {
413
414
415
416
        if (_rec_types.empty())
            return;

        for (size_t i = 0; i < _rec_types.size(); ++i)
417
418
419
420
        {
            if (_rec_types[i] != weight_type::REAL_NORMAL)
                continue;

421
422
423
424
            auto rec = _c_rec[i][e];
            auto d = (std::pow(rec, 2) -
                      std::pow(rec - delta[i], 2));
            _c_drec[i][e] += d;
425
426
427
        }
    }

428
429
    void remove_partition_node(size_t v, size_t r)
    {
430
431
        assert(size_t(_b[v]) == r);

432
433
434
435
436
437
438
439
440
441
442
443
444
        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);
            }
        }

445
        _wr[r] -= _vweight[v];
446

447
        if (!_egroups.empty() && _egroups_enabled)
448
            _egroups.remove_vertex(v, _b, _g);
449
450

        if (is_partition_stats_enabled())
451
452
453
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
454
455
456
457
    }

    void add_partition_node(size_t v, size_t r)
    {
458
459
        _b[v] = r;

460
        _wr[r] += _vweight[v];
461

462
        if (!_egroups.empty() && _egroups_enabled)
463
            _egroups.add_vertex(v, _b, _eweight, _g);
464
465

        if (is_partition_stats_enabled())
466
467
            get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
                                              _eweight, _degs);
Tiago Peixoto's avatar
Tiago Peixoto committed
468

469
        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
Tiago Peixoto's avatar
Tiago Peixoto committed
470
        {
471
            remove_element(_empty_blocks, _empty_pos, r);
472
            add_element(_candidate_blocks, _candidate_pos, r);
473
474
475
476
477
478
479

            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]);
            }
480
481
482
        }
    }

483
484
485
    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
486
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt));
487
488
489
490
    }

    void remove_vertex(size_t v, size_t r)
    {
491
        remove_vertex(v, r, [](auto&) { return false; });
492
493
    }

494
495
    void remove_vertex(size_t v)
    {
496
497
        size_t r = _b[v];
        remove_vertex(v, r);
498
499
    }

500
501
    template <class Vlist>
    void remove_vertices(Vlist& vs)
502
503
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
504
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
505
506
507

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

509
510
511
512
513
514
515
516
517
518
519
        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)
520
            remove_vertex(v, _b[v],
521
                          [&](auto& e) { return eset.find(e) != eset.end(); });
522
523
524
525
526
527
528
529

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

530
            auto me = _emat.get_me(r, s);
531
532
533
534
535
536
537
538
539

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

            assert(_mrs[me] >= 0);

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

540
            for (size_t i = 0; i < _rec_types.size(); ++i)
541
            {
542
543
544
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
545
                    _bdrec[i][me] -= _drec[i][e];
546
                    [[gnu::fallthrough]];
547
548
                default:
                    _brec[i][me] -= _rec[i][e];
549
                }
550
            }
551

552
553
554
555
556
            if (_mrs[me] == 0)
            {
                _emat.remove_me(me, _bg);
                if (_coupled_state != nullptr)
                    _coupled_state->remove_edge(me);
557
558
                else
                    boost::remove_edge(me, this->_bg);
559
            }
560
561
562
        }
    }

563
564
565
566
567
568
    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

569
570
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
571
    {
572
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
573
574
    }

575
576
    void add_vertex(size_t v, size_t r)
    {
577
        add_vertex(v, r, [](auto&){ return false; });
578
579
    }

580
581
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
582
    {
583
584
585
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
        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)
607
            add_vertex(vr.first, vr.second,
608
                       [&](auto& e){ return eset.find(e) != eset.end(); });
609
610
611
612
613
614
615
616
617
618

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

619
            if (me == _emat.get_null_edge())
620
            {
621
                me = boost::add_edge(r, s, _bg).first;
622
623
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
624
625
626
627
628
                for (size_t i = 0; i < _rec_types.size(); ++i)
                {
                    _c_brec[i][me] = 0;
                    _c_bdrec[i][me] = 0;
                }
629
630
631

                if (_coupled_state != nullptr)
                    _coupled_state->add_edge(me);
632
633
634
635
636
637
638
639
640
            }

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

            auto ew = _eweight[e];

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

642
            for (size_t i = 0; i < _rec_types.size(); ++i)
643
            {
644
645
646
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
647
                    _bdrec[i][me] += _drec[i][e];
648
                    [[gnu::fallthrough]];
649
650
                default:
                    _brec[i][me] += _rec[i][e];
651
                }
652
653
654
655
            }
        }
    }

656
    void add_vertices(python::object ovs, python::object ors)
657
    {
658
659
660
661
662
        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);
    }

663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
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
    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())++;
Tiago Peixoto's avatar
Tiago Peixoto committed
784
                if constexpr (is_directed_::apply<g_t>::type::value)
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
                    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())--;
Tiago Peixoto's avatar
Tiago Peixoto committed
801
                if constexpr (is_directed_::apply<g_t>::type::value)
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
                    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);
    }

821
822
823
824
825
826
827
828
829
830
831
832
833
    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)
    {
834
        _N -= vweight[v];
835
        vweight[v] = w;
836
        _N += w;
837
838
    }

839
840
841
842
843
844
845
846
847
848
849
850
851
    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));
852
        vweight[v] = 0;
853
854
    }

855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
    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);
    }

871
872
873
874
    template <class VMap>
    void set_partition(VMap&& b)
    {
        for (auto v : vertices_range(_g))
875
876
877
878
879
880
        {
            size_t r = b[v];
            while (r >= num_vertices(_bg))
                add_block();
            move_vertex(v, r);
        }
881
882
883
884
885
886
887
888
    }

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

889
890
891
892
893
894
895
896
897
898
899
900
901
902
    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>
903
    void merge_vertices(size_t u, size_t v, Emap&& ec)
904
    {
905
        merge_vertices(u, v, ec, _is_weighted);
906
907
    }

908
    template <class Emap>
909
    void merge_vertices(size_t, size_t, Emap&&, std::false_type)
910
911
912
913
914
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
915
    void merge_vertices(size_t u, size_t v, Emap&& ec, std::true_type)
916
    {
917
918
919
        if (u == v)
            return;

920
        auto eweight_c = _eweight.get_checked();
921
922
923
924
925
926
        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());
927
928
929
930
931
932

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

937
        size_t nrec = this->_rec_types.size();
938
        std::vector<double> ecc(nrec), decc(nrec);
939

940
941
942
943
944
945
946
        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;
947

948
949
            std::fill(ecc.begin(), ecc.end(), 0);
            std::fill(decc.begin(), decc.end(), 0);
950
            for (auto& e : es)
951
            {
952
                w += _eweight[e];
953
954
955
956
957
                for (size_t i = 0; i < nrec; ++i)
                {
                    ecc[i] += _rec[i][e];
                    decc[i] += _drec[i][e];
                }
958
            }
959
960
961
962

            if (t == u)
            {
                t = v;
Tiago Peixoto's avatar
Tiago Peixoto committed
963
                if constexpr (!is_directed_::apply<g_t>::type::value)
964
965
966
                {
                    assert(w % 2 == 0);
                    w /= 2;
967
968
969
970
971
                    for (size_t i = 0; i < nrec; ++i)
                    {
                        ecc[i] /= 2;
                        decc[i] /= 2;
                    }
972
973
974
                }
            }

975
            auto iter = ns_v.find(std::make_tuple(t, l));
976
977
978
979
            if (iter != ns_v.end())
            {
                auto& e = iter->second.front();
                _eweight[e] += w;
980
981
982
983
984
                for (size_t i = 0; i < nrec; ++i)
                {
                    _rec[i][e] += ecc[i];
                    _drec[i][e] += decc[i];
                }
985
                assert(ec[e] == l);
986
987
988
            }
            else
            {
989
                auto e = boost::add_edge(v, t, _g).first;
990
                ns_v[std::make_tuple(t, l)].push_back(e);
991
                eweight_c[e] = w;
992
993
994
995
996
                for (size_t i = 0; i < nrec; ++i)
                {
                    c_rec[i][e] = ecc[i];
                    c_drec[i][e] = decc[i];
                }
997
                set_prop(ec, e, l);
998
                assert(ec[e] == l);
999
1000
            }
        }
For faster browsing, not all history is shown. View entire blame