graph_blockmodel.hh 90.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-2017 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
#define BLOCK_STATE_params                                                     \
    ((g, &, all_graph_views, 1))                                               \
63
64
    ((is_weighted,, bool_tr, 1))                                               \
    ((use_hash,, bool_tr, 1))                                                  \
65
    ((_abg, &, boost::any&, 0))                                                \
66
67
    ((_aeweight, &, boost::any&, 0))                                           \
    ((_avweight, &, boost::any&, 0))                                           \
68
    ((_adegs, &, boost::any&, 0))                                              \
69
70
71
72
73
    ((mrs,, emap_t, 0))                                                        \
    ((mrp,, vmap_t, 0))                                                        \
    ((mrm,, vmap_t, 0))                                                        \
    ((wr,, vmap_t, 0))                                                         \
    ((b,, vmap_t, 0))                                                          \
74
75
76
77
    ((empty_blocks, & ,std::vector<size_t>&, 0))                               \
    ((empty_pos,, vmap_t, 0))                                                  \
    ((candidate_blocks, &, std::vector<size_t>&, 0))                           \
    ((candidate_pos,, vmap_t, 0))                                              \
78
79
80
81
    ((bclabel,, vmap_t, 0))                                                    \
    ((pclabel,, vmap_t, 0))                                                    \
    ((merge_map,, vmap_t, 0))                                                  \
    ((deg_corr,, bool, 0))                                                     \
82
    ((rec_types,, std::vector<int32_t>, 0))                                    \
83
84
85
86
    ((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))                      \
87
    ((brecsum,, vprop_map_t<double>::type, 0))                                 \
88
    ((wparams, &, std::vector<std::vector<double>>&, 0))                       \
89
90
    ((recdx, &, std::vector<double>&, 0))                                      \
    ((Lrecdx, &, std::vector<double>&, 0))                                     \
91
    ((epsilon, &, std::vector<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

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

template <class... Ts>
class BlockState
100
    : public BlockStateBase<Ts...>, public BlockStateVirtualBase
101
102
103
104
105
106
107
108
109
110
111
{
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
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
114
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
115
          _emat(_bg, rng),
116
          _neighbor_sampler(_g, _eweight),
117
          _m_entries(num_vertices(_bg))
118
    {
119
120
        _empty_blocks.clear();
        _candidate_blocks.clear();
121
        _candidate_blocks.push_back(null_group);
122
123
124
125
126
127
128
        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);
        }
129
130
131
132
        for (auto& p : _rec)
            _c_rec.push_back(p.get_checked());
        for (auto& p : _drec)
            _c_drec.push_back(p.get_checked());
133
134
135
136
137
138
139
140
141
142
143
        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());

144
        if (!_rec_types.empty())
145
        {
146
147
            _recx2.resize(_rec_types.size());
            _recdx.resize(_rec_types.size());
148
            for (auto me : edges_range(_bg))
149
            {
150
151
152
153
154
                if (_brec[0][me] > 0)
                {
                    _B_E++;
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                    {
155
                        if (_rec_types[i] == weight_type::REAL_NORMAL)
156
157
158
159
160
161
162
163
164
165
166
                        {
                            _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++;
167
168
            }
        }
169
170
171
172
173
174
175
176
177

        _rt = weight_type::NONE;
        for (auto rt : _rec_types)
        {
            _rt = rt;
            if (rt == weight_type::REAL_NORMAL)
                break;
        }
        _dBdx.resize(_rec_types.size());
178
179
180
181
182
183
    }

    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()),
184
185
          _c_rec(other._c_rec),
          _c_drec(other._c_drec),
186
187
188
189
          _c_brec(other._c_brec),
          _c_bdrec(other._c_bdrec),
          _recsum(other._recsum),
          _recx2(other._recx2),
190
          _dBdx(other._dBdx),
191
          _B_E(other._B_E),
192
193
          _B_E_D(other._B_E_D),
          _rt(other._rt),
194
195
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
196
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
197
          _emat(other._emat),
198
          _egroups_enabled(other._egroups_enabled),
199
          _neighbor_sampler(other._neighbor_sampler),
200
          _m_entries(num_vertices(_bg))
201
202
203
204
205
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

206
207
208
    // =========================================================================
    // State modification
    // =========================================================================
209

210
    template <class MEntries, class EFilt>
211
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
212
                          EFilt&& efilt)
213
214
    {
        auto mv_entries = [&](auto&&... args)
215
            {
216
217
218
                move_entries(v, r, nr, _b, _g, _eweight, num_vertices(_bg),
                             m_entries, std::forward<EFilt>(efilt),
                             is_loop_nop(),
219
                             std::forward<decltype(args)>(args)...);
220
            };
221

222
        if (_rt == weight_type::NONE)
223
        {
224
            mv_entries();
225
        }
226
        else
227
        {
228
229
230
231
            if (_rt == weight_type::REAL_NORMAL)
                mv_entries(_rec, _drec);
            else
                mv_entries(_rec);
232
233
        }
    }
234

235
236
237
    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
238
        get_move_entries(v, r, nr, m_entries, [](auto) {return false;});
239
    }
240
241


242
243
    template <bool Add, class EFilt>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt)
244
245
    {
        if (Add)
246
            get_move_entries(v, null_group, r, _m_entries,
247
                             std::forward<EFilt>(efilt));
248
        else
249
            get_move_entries(v, r, null_group, _m_entries,
250
                             std::forward<EFilt>(efilt));
Tiago Peixoto's avatar
Tiago Peixoto committed
251

252
253
254
255
        auto eops = [&](auto&& mid_op, auto&& end_op)
            {
                entries_op(_m_entries, _emat,
                           [&](auto r, auto s, auto& me, auto& delta)
256
                           {
257
258
                               if (get<0>(delta) == 0) // can happen with
                                   return;             // zero-weight edges
259

260
261
262
                               if (Add && me == _emat.get_null_edge())
                               {
                                   me = boost::add_edge(r, s, this->_bg).first;
263
                                   this->_emat.put_me(r, s, me);
264
265
266
267
268
269
270
                                   this->_c_mrs[me] = 0;
                                   for (size_t i = 0; i < this->_rec_types.size(); ++i)
                                   {
                                       this->_c_brec[i][me] = 0;
                                       this->_c_bdrec[i][me] = 0;
                                   }
                               }
271

272
                               mid_op(me, delta);
273

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

278
279
280
                               assert(this->_mrs[me] >= 0);
                               assert(this->_mrp[r] >= 0);
                               assert(this->_mrm[s] >= 0);
281

282
283
284
285
286
287
                               end_op(me, delta);
                           });
                   };

        if (_rec_types.empty())
        {
288
            eops([](auto&, auto&){}, [](auto&, auto&){});
289
290
291
292
293
294
        }
        else
        {
            auto end_op = [&](auto& me, auto& delta)
                {
                    for (size_t i = 0; i < this->_rec_types.size(); ++i)
295
296
297
298
299
300
301
302
303
304
                    {
                        switch (this->_rec_types[i])
                        {
                        case weight_type::REAL_NORMAL: // signed weights
                            this->_bdrec[i][me] += get<2>(delta)[i];
                            [[gnu::fallthrough]];
                        default:
                            this->_brec[i][me] += get<1>(delta)[i];
                        }
                    }
305
                };
306

307
308
309
            auto mid_op_BE =
                [&](auto& me, auto&& delta)
                {
310
                    auto mrs = this->_brec[0][me];
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
                    if (Add && mrs == 0 && mrs + get<1>(delta)[0] > 0)
                    {
                        _B_E++;
                        if (_coupled_state != nullptr)
                            _coupled_state->add_edge(me);
                    }

                    if (!Add && mrs > 0 && mrs + get<1>(delta)[0] == 0)
                    {
                        _B_E--;
                        if (_coupled_state != nullptr)
                            _coupled_state->remove_edge(me);
                    }
                };

            if (_rt != weight_type::REAL_NORMAL)
            {
                eops(mid_op_BE, end_op);
            }
            else
            {
                auto mid_op =
                    [&](auto& me, auto&& delta)
                    {
                        auto& mrs = this->_brec[0][me];
                        mid_op_BE(me, delta);

338
339
340
                        auto n_mrs = mrs + get<1>(delta)[0];

                        if (n_mrs > 1)
341
                        {
342
                            if (Add && mrs < 2)
343
344
345
                            {
                                if (_B_E_D == 0 && this->_Lrecdx[0] >= 0)
                                    this->_Lrecdx[0] += 1;
346
                                _B_E_D++;
347
                            }
348
349
350

                            for (size_t i = 0; i < this->_rec_types.size(); ++i)
                            {
351
352
353
354
355
356
357
                                if (this->_rec_types[i] != weight_type::REAL_NORMAL)
                                    continue;
                                auto dx = \
                                    (this->_bdrec[i][me] + get<2>(delta)[i]
                                     - (std::pow((this->_brec[i][me] +
                                                  get<1>(delta)[i]), 2) / n_mrs));
                                this->_recdx[i] += dx;
358
                            }
359
360
361
362
363
                        }

                        if (mrs > 1)
                        {
                            if (!Add && n_mrs < 2)
364
                            {
365
                                _B_E_D--;
366
367
368
                                if (_B_E_D == 0 && this->_Lrecdx[0] >= 0)
                                    this->_Lrecdx[0] -= 1;
                            }
369
370
371
372
373
374
375
376
377

                            for (size_t i = 0; i < this->_rec_types.size(); ++i)
                            {
                                if (this->_rec_types[i] != weight_type::REAL_NORMAL)
                                    continue;
                                auto dx = (this->_bdrec[i][me] -
                                           std::pow(this->_brec[i][me], 2) / mrs);
                                this->_recdx[i] -= dx;
                            }
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
                        }

                        for (size_t i = 0; i < this->_rec_types.size(); ++i)
                        {
                            if (this->_rec_types[i] == weight_type::REAL_NORMAL)
                            {
                                _recx2[i] -= std::pow(this->_brec[i][me], 2);
                                _recx2[i] += std::pow(this->_brec[i][me] +
                                                      get<1>(delta)[i], 2);
                            }
                        }

                    };

                auto coupled_end_op = [&](auto& me, auto& delta)
                    {
                        end_op(me, delta);
                        if (_coupled_state != nullptr)
                            _coupled_state->update_edge(me, get<1>(delta));
                    };

399
400
                if (_Lrecdx[0] >= 0)
                {
401
402
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                        _Lrecdx[i+1] -= _recdx[i] * _B_E_D;
403
                }
404

405
                eops(mid_op, coupled_end_op);
406

407
408
                if (_Lrecdx[0] >= 0)
                {
409
410
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                        _Lrecdx[i+1] += _recdx[i] * _B_E_D;
411
                }
412
413
            }
        }
414

415
        if (!_rec_types.empty() &&
416
            _rec_types[1] == weight_type::DELTA_T) // waiting times
417
        {
418
            if (_ignore_degrees[v] > 0)
419
            {
420
                auto dt = out_degreeS()(v, _g, _rec[1]);
421
                if (Add)
422
                    _brecsum[r] += dt;
423
                else
424
                    _brecsum[r] -= dt;
425
            }
426
427
        }

428
        if (Add)
429
            BlockState::add_partition_node(v, r);
430
        else
431
            BlockState::remove_partition_node(v, r);
432
433
    }

434
435
    void add_edge(const GraphInterface::edge_t& e)
    {
436
437
438
439
440
        if (_rec_types.empty())
            return;
        auto crec = _rec[0].get_checked();
        crec[e] = 1;
        for (size_t i = 1; i < _rec_types.size(); ++i)
441
        {
442
443
444
445
446
447
448
449
            auto drec = _drec[i].get_checked();
            drec[e] = 0;
        }
        size_t r = _b[source(e, _g)];
        size_t s = _b[target(e, _g)];
        auto& me = _emat.get_me(r, s);
        assert(me != _emat.get_null_edge());
        _brec[0][me] += 1;
450

451
452
        if (_brec[0][me] == 1)
            _B_E++;
453

454
455
456
457
458
459
460
        size_t old_B_E_D = _B_E_D;
        if (_brec[0][me] == 2)
        {
            if (_B_E_D == 0 && _Lrecdx[0] >= 0)
                _Lrecdx[0] += 1;
            _B_E_D++;
        }
461

462
463
464
        if (_rt == weight_type::REAL_NORMAL && _brec[0][me] > 1)
        {
            for (size_t i = 0; i < _rec_types.size(); ++i)
465
            {
466
467
                if (_rec_types[i] != weight_type::REAL_NORMAL)
                    continue;
468

469
470
                if (_Lrecdx[0] >= 0)
                    _Lrecdx[i+1] -= _recdx[i] * old_B_E_D;
471

472
473
474
475
476
477
478
479
                auto dx = (_bdrec[i][me] -
                           std::pow(_brec[i][me], 2) / _brec[0][me]);
                _recdx[i] += dx;
                if (_brec[0][me] > 2)
                {
                    dx = (_bdrec[i][me] -
                          std::pow(_brec[i][me], 2) / (_brec[0][me] - 1));
                    _recdx[i] -= dx;
480
                }
481
482
483

                if (_Lrecdx[0] >= 0)
                    _Lrecdx[i+1] += _recdx[i] * _B_E_D;
484
            }
485
486
487
488
489
        }
    }

    void remove_edge(const GraphInterface::edge_t& e)
    {
490
491
        if (_rec_types.empty())
            return;
492

493
494
495
496
497
        size_t r = _b[source(e, _g)];
        size_t s = _b[target(e, _g)];
        auto& me = _emat.get_me(r, s);
        _brec[0][me] -= 1;
        _rec[0][e] = 0;
498

499
500
501
502
503
504
505
506
507
508
        if (_brec[0][me] == 0)
            _B_E--;

        size_t old_B_E_D = _B_E_D;
        if (_brec[0][me] == 1)
        {
            _B_E_D--;
            if (_B_E_D == 0 && _Lrecdx[0] >= 0)
                _Lrecdx[0] -= 1;
        }
509

510
511
512
        if (_rt == weight_type::REAL_NORMAL && _brec[0][me] > 0)
        {
            for (size_t i = 0; i < _rec_types.size(); ++i)
513
            {
514
515
516
517
518
519
520
521
                if (_rec_types[i] != weight_type::REAL_NORMAL)
                    continue;
                if (_Lrecdx[0] >= 0)
                    _Lrecdx[i+1] -= _recdx[i] * old_B_E_D;
                auto dx = (_bdrec[i][me] -
                           std::pow(_brec[i][me], 2) / (_brec[0][me] + 1));
                _recdx[i] -= dx;
                if (_brec[0][me] > 1)
522
                {
523
524
525
                    dx = (_bdrec[i][me] -
                          std::pow(_brec[i][me], 2) / _brec[0][me]);
                    _recdx[i] += dx;
526
                }
527
528
                if (_Lrecdx[0] >= 0)
                    _Lrecdx[i+1] += _recdx[i] * _B_E_D;
529
            }
530
531
532
        }
    }

533
534
535
    void update_edge(const GraphInterface::edge_t& e,
                     const std::vector<double>& delta)
    {
536
537
538
        if (_rec_types.empty())
            return;

539
540
541
        size_t r = _b[source(e, _g)];
        size_t s = _b[target(e, _g)];
        auto& me = _emat.get_me(r, s);
542
        auto ers = _brec[0][me];
543
        for (size_t i = 0; i < _rec_types.size(); ++i)
544
545
546
547
        {
            if (_rec_types[i] != weight_type::REAL_NORMAL)
                continue;

548
549
550
551
            auto rec = _c_rec[i][e];
            auto d = (std::pow(rec, 2) -
                      std::pow(rec - delta[i], 2));
            _c_drec[i][e] += d;
552
553
554
555
            _bdrec[i][me] += d;
            if (ers > 1)
            {
                _recdx[i] += d;
556
557
                if (_Lrecdx[0] >= 0)
                    _Lrecdx[i+1] += d * _B_E_D;
558
559
560
561
            }
        }
    }

562
563
    void remove_partition_node(size_t v, size_t r)
    {
564
565
        assert(size_t(_b[v]) == r);

566
        _wr[r] -= _vweight[v];
567

568
        if (!_egroups.empty() && _egroups_enabled)
569
            _egroups.remove_vertex(v, _b, _g);
570
571

        if (is_partition_stats_enabled())
572
573
574
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
575

576
        if (_vweight[v] > 0 && _wr[r] == 0)
577
        {
578
            remove_element(_candidate_blocks, _candidate_pos, r);
579
580
581
582
583
584
            add_element(_empty_blocks, _empty_pos, r);
        }
    }

    void add_partition_node(size_t v, size_t r)
    {
585
586
        _b[v] = r;

587
        _wr[r] += _vweight[v];
588

589
        if (!_egroups.empty() && _egroups_enabled)
590
            _egroups.add_vertex(v, _b, _eweight, _g);
591
592

        if (is_partition_stats_enabled())
593
594
            get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
                                              _eweight, _degs);
Tiago Peixoto's avatar
Tiago Peixoto committed
595

596
        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
Tiago Peixoto's avatar
Tiago Peixoto committed
597
        {
598
            remove_element(_empty_blocks, _empty_pos, r);
599
            add_element(_candidate_blocks, _candidate_pos, r);
600
601
602
        }
    }

603
604
605
    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
606
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt));
607
608
609
610
611
    }

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

614
615
    void remove_vertex(size_t v)
    {
616
617
        size_t r = _b[v];
        remove_vertex(v, r);
618
619
    }

620
621
    template <class Vlist>
    void remove_vertices(Vlist& vs)
622
623
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
624
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
625
626
627

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

629
630
631
632
633
634
635
636
637
638
639
        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)
640
            remove_vertex(v, _b[v],
641
                          [&](auto& e) { return eset.find(e) != eset.end(); });
642
643
644
645
646
647
648
649

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

650
            auto me = _emat.get_me(r, s);
651
652
653
654
655
656
657
658
659

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

            assert(_mrs[me] >= 0);

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

660
            for (size_t i = 0; i < _rec_types.size(); ++i)
661
            {
662
663
664
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
665
                    _bdrec[i][me] -= _drec[i][e];
666
                    [[gnu::fallthrough]];
667
668
                default:
                    _brec[i][me] -= _rec[i][e];
669
                }
670
            }
671

672
673
            // if (_mrs[me] == 0)
            //     _emat.remove_me(me, _bg);
674
675
676
        }
    }

677
678
679
680
681
682
    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

683
684
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
685
    {
686
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
687
688
    }

689
690
    void add_vertex(size_t v, size_t r)
    {
691
        add_vertex(v, r, [](auto&){ return false; });
692
693
    }

694
695
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
696
    {
697
698
699
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
        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)
721
            add_vertex(vr.first, vr.second,
722
                       [&](auto& e){ return eset.find(e) != eset.end(); });
723
724
725
726
727
728
729
730
731
732

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

733
            if (me == _emat.get_null_edge())
734
            {
735
                me = boost::add_edge(r, s, _bg).first;
736
737
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
738
739
740
741
742
                for (size_t i = 0; i < _rec_types.size(); ++i)
                {
                    _c_brec[i][me] = 0;
                    _c_bdrec[i][me] = 0;
                }
743
744
745
746
747
748
749
750
751
            }

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

            auto ew = _eweight[e];

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

753
            for (size_t i = 0; i < _rec_types.size(); ++i)
754
            {
755
756
757
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
758
                    _bdrec[i][me] += _drec[i][e];
759
                    [[gnu::fallthrough]];
760
761
                default:
                    _brec[i][me] += _rec[i][e];
762
                }
763
764
765
766
            }
        }
    }

767
    void add_vertices(python::object ovs, python::object ors)
768
    {
769
770
771
772
773
774
775
776
777
778
        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
779
            return _bclabel[r] == _bclabel[nr];
780
781
    }

782
    // move a vertex from its current block to block nr
783
    void move_vertex(size_t v, size_t r, size_t nr)
784
785
786
    {
        if (r == nr)
            return;
787
788

        if (!allow_move(r, nr))
789
            throw ValueException("cannot move vertex across clabel barriers");
790
791
792

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

799
            if (_wr[nr] == 0)
800
801
802
803
804
805
            {
                _coupled_state->set_vertex_weight(nr, 1);
                _coupled_state->add_partition_node(nr, _bclabel[r]);
                _bclabel[nr] = _bclabel[r];
            }
        }
806
807
808

        remove_vertex(v, r, [](auto&) {return false;});
        add_vertex(v, nr, [](auto&) {return false;});
Tiago Peixoto's avatar
Tiago Peixoto committed
809
810
    }

811
812
    void move_vertex(size_t v, size_t nr)
    {
813
814
        size_t r = _b[v];
        move_vertex(v, r, nr);
815
816
    }

817
818
819
820
821
822
823
824
825
826
827
828
829
830
    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;
831
832
    }

833
834
835
836
837
838
839
840
841
842
843
844
845
    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));
846
        vweight[v] = 0;
847
848
    }

849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
    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);
    }

865
866
867
868
869
870
871
872
873
874
875
876
877
    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());
    }

878
879
880
881
882
883
884
885
886
887
888
889
890
891
    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>
892
    void merge_vertices(size_t u, size_t v, Emap&& ec)
893
    {
894
        merge_vertices(u, v, ec, _is_weighted);
895
896
    }

897
    template <class Emap>
898
    void merge_vertices(size_t, size_t, Emap&&, std::false_type)
899
900
901
902
903
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
904
    void merge_vertices(size_t u, size_t v, Emap&& ec, std::true_type)
905
    {
906
907
908
        if (u == v)
            return;

909
        auto eweight_c = _eweight.get_checked();
910
911
912
913
914
915
        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());
916
917
918
919
920
921

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

926
        size_t nrec = this->_rec_types.size();
927
        std::vector<double> ecc(nrec), decc(nrec);
928

929
930
931
932
933
934
935
        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;
936

937
938
            std::fill(ecc.begin(), ecc.end(), 0);
            std::fill(decc.begin(), decc.end(), 0);
939
            for (auto& e : es)
940
            {
941
                w += _eweight[e];
942
943
944
945
946
                for (size_t i = 0; i < nrec; ++i)
                {
                    ecc[i] += _rec[i][e];
                    decc[i] += _drec[i][e];
                }
947
            }
948
949
950
951
952
953
954
955

            if (t == u)
            {
                t = v;
                if (!is_directed::apply<g_t>::type::value)
                {
                    assert(w % 2 == 0);
                    w /= 2;
956
957
958
959
960
                    for (size_t i = 0; i < nrec; ++i)
                    {
                        ecc[i] /= 2;
                        decc[i] /= 2;
                    }
961
962
963
                }
            }

964
            auto iter = ns_v.find(std::make_tuple(t, l));
965
966
967
968
            if (iter != ns_v.end())
            {
                auto& e = iter->second.front();
                _eweight[e] += w;
969
970
971
972
973
                for (size_t i = 0; i < nrec; ++i)
                {
                    _rec[i][e] += ecc[i];
                    _drec[i][e] += decc[i];
                }
974
                assert(ec[e] == l);
975
976
977
            }
            else
            {
978
                auto e = boost::add_edge(v, t, _g).first;
979
                ns_v[std::make_tuple(t, l)].push_back(e);
980
                eweight_c[e] = w;
981
982
983
984
985
                for (size_t i = 0; i < nrec; ++i)
                {
                    c_rec[i][e] = ecc[i];
                    c_drec[i][e] = decc[i];
                }
986
                set_prop(ec, e, l);
987
                assert(ec[e] == l);
988
989
990
991
992
993
994
995
996
            }
        }

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

            for(auto e : in_edges_range(v, _g))
997
                ns_v[std::make_tuple(source(e, _g), ec[e])].push_back(e);
998
            for(auto e : in_edges_range(u, _g))
999
                ns_u[std::make_tuple(source(e, _g), ec[e])].push_back(e);
1000

For faster browsing, not all history is shown. View entire blame