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

#ifndef GRAPH_BLOCKMODEL_HH
#define GRAPH_BLOCKMODEL_HH

#include "config.h"

#include <vector>

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

namespace graph_tool
{
using namespace boost;
using namespace std;

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

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

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

enum weight_type
{
    NONE,
51 52
    REAL_EXPONENTIAL,
    REAL_NORMAL,
53 54
    DISCRETE_GEOMETRIC,
    DISCRETE_POISSON,
55
    DISCRETE_BINOMIAL,
56 57 58
    DELTA_T
};

59 60 61 62

#define BLOCK_STATE_params                                                     \
    ((g, &, all_graph_views, 1))                                               \
    ((degs,, degs_tr, 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 69 70 71 72
    ((mrs,, emap_t, 0))                                                        \
    ((mrp,, vmap_t, 0))                                                        \
    ((mrm,, vmap_t, 0))                                                        \
    ((wr,, vmap_t, 0))                                                         \
    ((b,, vmap_t, 0))                                                          \
73 74 75 76
    ((empty_blocks, & ,std::vector<size_t>&, 0))                               \
    ((empty_pos,, vmap_t, 0))                                                  \
    ((candidate_blocks, &, std::vector<size_t>&, 0))                           \
    ((candidate_pos,, vmap_t, 0))                                              \
77 78 79 80
    ((bclabel,, vmap_t, 0))                                                    \
    ((pclabel,, vmap_t, 0))                                                    \
    ((merge_map,, vmap_t, 0))                                                  \
    ((deg_corr,, bool, 0))                                                     \
81 82 83 84 85
    ((rec_types,, std::vector<int32_t>, 0))                                    \
    ((rec,, eprop_map_t<std::vector<double>>::type, 0))                        \
    ((drec,, eprop_map_t<std::vector<double>>::type, 0))                       \
    ((brec,, eprop_map_t<std::vector<double>>::type, 0))                       \
    ((bdrec,, eprop_map_t<std::vector<double>>::type, 0))                      \
86
    ((brecsum,, vprop_map_t<double>::type, 0))                                 \
87
    ((wparams, &, std::vector<std::vector<double>>&, 0))                       \
88
    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))                \
89
    ((bignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))               \
90
    ((allow_empty,, bool, 0))
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

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

    template <class RNG, class... ATs,
              typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
    BlockState(RNG& rng, ATs&&... args)
        : BlockStateBase<Ts...>(std::forward<ATs>(args)...),
          _bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
          _c_mrs(_mrs.get_checked()),
108 109 110 111
          _c_brec(_brec.get_checked()),
          _c_bdrec(_bdrec.get_checked()),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
112
          _emat(_bg, rng),
113
          _egroups_enabled(true),
114
          _neighbour_sampler(_g, _eweight),
115
          _m_entries(num_vertices(_bg)),
116
          _coupled_state(nullptr)
117
    {
118 119
        _empty_blocks.clear();
        _candidate_blocks.clear();
120
        _candidate_blocks.push_back(null_group);
121 122 123 124 125 126 127
        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);
        }
128 129 130 131 132 133
    }

    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()),
134 135 136 137
          _c_brec(_brec.get_checked()),
          _c_bdrec(_bdrec.get_checked()),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
138
          _emat(other._emat),
139
          _egroups_enabled(other._egroups_enabled),
140
          _neighbour_sampler(other._neighbour_sampler),
141
          _m_entries(num_vertices(_bg)),
142
          _coupled_state(nullptr)
143 144 145 146 147
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

148 149 150
    // =========================================================================
    // State modification
    // =========================================================================
151

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

163 164
        int rec_type = weight_type::NONE;
        for (auto rt : _rec_types)
165
        {
166 167 168 169 170 171 172 173
            rec_type = rt;
            if (rt == weight_type::REAL_NORMAL)
                break;
        }

        switch (rec_type)
        {
        case weight_type::REAL_EXPONENTIAL:
174 175
        case weight_type::DISCRETE_GEOMETRIC:
        case weight_type::DISCRETE_POISSON:
176
        case weight_type::DISCRETE_BINOMIAL:
177
        case weight_type::DELTA_T:
178
            mv_entries(_rec);
179
            break;
180
        case weight_type::REAL_NORMAL:
181
            mv_entries(_rec, _drec);
182 183 184 185 186
            break;
        default: // no weights
            mv_entries();
        }
    }
187

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


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

206 207 208
        entries_op(_m_entries, _emat,
                   [&](auto r, auto s, auto& me, auto& delta)
                   {
209 210 211
                       if (get<0>(delta) == 0) // can happen with zero-weight
                           return;             // edges

212
                       if (Add && me == _emat.get_null_edge())
213 214 215 216
                       {
                           me = add_edge(r, s, this->_bg).first;
                           _emat.put_me(r, s, me);
                           this->_c_mrs[me] = 0;
217 218
                           this->_c_brec[me].clear();
                           this->_c_bdrec[me].clear();
219 220 221 222 223 224
                       }

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

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

229 230 231
                       this->_brec[me].resize(get<1>(delta).size());
                       this->_bdrec[me].resize(get<2>(delta).size());
                       for (size_t i = 0; i < this->_rec_types.size(); ++i)
232
                       {
233 234 235 236 237 238 239 240 241 242 243
                           switch (this->_rec_types[i])
                           {
                           case weight_type::REAL_NORMAL: // signed weights
                               this->_bdrec[me][i] += get<2>(delta)[i];
                           case weight_type::REAL_EXPONENTIAL:
                           case weight_type::DISCRETE_GEOMETRIC:
                           case weight_type::DISCRETE_POISSON:
                           case weight_type::DISCRETE_BINOMIAL:
                           case weight_type::DELTA_T:
                               this->_brec[me][i] += get<1>(delta)[i];
                           }
244
                       }
245 246 247 248 249

                       if (!Add && this->_mrs[me] == 0)
                       {
                           _emat.remove_me(me, this->_bg);
                       }
250 251
                   });

252 253
        if (!_rec_types.empty() &&
            _rec_types.front() == weight_type::DELTA_T) // waiting times
254
        {
255
            if (_ignore_degrees[v] > 0)
256
            {
257
                auto dt = out_degreeS()(v, _g, _rec);
258
                if (Add)
259
                    _brecsum[r] += dt[0];
260
                else
261
                    _brecsum[r] -= dt[0];
262
            }
263 264
        }

265
        if (Add)
266
        {
267
            _b[v] = r;
268
            add_partition_node(v, r);
269
        }
270
        else
271
        {
272
            remove_partition_node(v, r);
273
        }
274 275 276 277
    }

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

280
        if (!_egroups.empty() && _egroups_enabled)
281
            _egroups.remove_vertex(v, _b, _g);
282 283

        if (is_partition_stats_enabled())
284 285 286
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
287

288
        if (_vweight[v] > 0 && _wr[r] == 0)
289
        {
290
            remove_element(_candidate_blocks, _candidate_pos, r);
291 292 293 294 295 296
            add_element(_empty_blocks, _empty_pos, r);
        }
    }

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

299
        if (!_egroups.empty() && _egroups_enabled)
300
            _egroups.add_vertex(v, _b, _eweight, _g);
301 302

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

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

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

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

324 325
    void remove_vertex(size_t v)
    {
326 327
        size_t r = _b[v];
        remove_vertex(v, r);
328 329
    }

330 331
    template <class Vlist>
    void remove_vertices(Vlist& vs)
332 333
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
334
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349

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

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

        for (auto v : vset)
350
            remove_vertex(v, _b[v],
351
                          [&](auto& e) { return eset.find(e) != eset.end(); });
352 353 354 355 356 357 358 359

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

360
            auto me = _emat.get_me(r, s);
361 362 363 364 365 366 367 368 369

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

            assert(_mrs[me] >= 0);

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

370
            for (size_t i = 0; i < _rec_types.size(); ++i)
371
            {
372 373 374 375 376 377 378 379 380 381
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
                    _bdrec[me][i] -= _drec[e][i];
                case weight_type::REAL_EXPONENTIAL:
                case weight_type::DISCRETE_GEOMETRIC:
                case weight_type::DISCRETE_POISSON:
                case weight_type::DISCRETE_BINOMIAL:
                    _brec[me][i] -= _rec[e][i];
                }
382
            }
383 384 385

            if (_mrs[me] == 0)
                _emat.remove_me(me, _bg);
386 387 388
        }
    }

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

395 396
    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
397
    {
398
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
399 400
    }

401 402
    void add_vertex(size_t v, size_t r)
    {
403
        add_vertex(v, r, [](auto&){ return false; });
404 405
    }

406 407
    template <class Vlist, class Blist>
    void add_vertices(Vlist& vs, Blist& rs)
408
    {
409 410 411
        if (vs.size() != rs.size())
            throw ValueException("vertex and group lists do not have the same size");

412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432
        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)
433
            add_vertex(vr.first, vr.second,
434
                       [&](auto& e){ return eset.find(e) != eset.end(); });
435 436 437 438 439 440 441 442 443 444

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

445
            if (me == _emat.get_null_edge())
446 447 448 449
            {
                me = add_edge(r, s, _bg).first;
                _emat.put_me(r, s, me);
                _c_mrs[me] = 0;
450 451
                _c_brec[me].clear();
                _c_bdrec[me].clear();
452 453 454 455 456 457 458 459 460
            }

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

            auto ew = _eweight[e];

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

462
            for (size_t i = 0; i < _rec_types.size(); ++i)
463
            {
464 465 466 467 468 469 470 471 472 473
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
                    _bdrec[me][i] += _drec[e][i];
                case weight_type::REAL_EXPONENTIAL:
                case weight_type::DISCRETE_GEOMETRIC:
                case weight_type::DISCRETE_POISSON:
                case weight_type::DISCRETE_BINOMIAL:
                    _brec[me][i] += _rec[e][i];
                }
474 475 476 477
            }
        }
    }

478
    void add_vertices(python::object ovs, python::object ors)
479
    {
480 481 482 483 484 485 486 487 488 489
        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
490
            return _bclabel[r] == _bclabel[nr];
491 492
    }

493
    // move a vertex from its current block to block nr
494
    void move_vertex(size_t v, size_t r, size_t nr)
495 496 497
    {
        if (r == nr)
            return;
498 499

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

502 503
        remove_vertex(v, r, [](auto&) {return false;});
        add_vertex(v, nr, [](auto&) {return false;});
504 505 506 507 508 509 510 511 512 513 514 515 516

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

            if (_wr[nr] == _vweight[v])
            {
                _coupled_state->set_vertex_weight(nr, 1);
                _coupled_state->add_partition_node(nr, _bclabel[r]);
517
                _coupled_state->_b[nr] = _bclabel[r];
518 519 520
                _bclabel[nr] = _bclabel[r];
            }
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
521 522
    }

523 524
    void move_vertex(size_t v, size_t nr)
    {
525 526
        size_t r = _b[v];
        move_vertex(v, r, nr);
527 528
    }

529 530 531 532 533 534 535 536 537 538 539 540 541 542
    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;
543 544
    }

545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
    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);
    }

561 562 563 564 565 566 567 568 569 570 571 572 573
    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());
    }

574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
    size_t virtual_remove_size(size_t v)
    {
        return _wr[_b[v]] - _vweight[v];
    }

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

    template <class Emap>
    void merge_vertices(size_t u, size_t v, Emap& ec)
    {
        if (u == v)
            return;
592
        merge_vertices(u, v, ec, _is_weighted);
593 594
    }

595 596
    template <class Emap>
    void merge_vertices(size_t, size_t, Emap&, std::false_type)
597 598 599 600 601
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

    template <class Emap>
602
    void merge_vertices(size_t u, size_t v, Emap& ec, std::true_type)
603 604
    {
        auto eweight_c = _eweight.get_checked();
605 606
        auto rec_c = _rec.get_checked();
        auto drec_c = _drec.get_checked();
607 608 609 610 611 612

        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))
613
            ns_u[std::make_tuple(target(e, _g), ec[e])].push_back(e);
614
        for(auto e : out_edges_range(v, _g))
615
            ns_v[std::make_tuple(target(e, _g), ec[e])].push_back(e);
616 617 618 619 620 621 622 623

        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;
624
            std::vector<double> ecc, decc;
625
            for (auto& e : es)
626
            {
627
                w += _eweight[e];
628 629 630
                ecc += _rec[e];
                decc += _drec[e];
            }
631 632 633 634 635 636 637 638

            if (t == u)
            {
                t = v;
                if (!is_directed::apply<g_t>::type::value)
                {
                    assert(w % 2 == 0);
                    w /= 2;
639 640
                    ecc /= 2;
                    decc /= 2;
641 642 643
                }
            }

644
            auto iter = ns_v.find(std::make_tuple(t, l));
645 646 647 648
            if (iter != ns_v.end())
            {
                auto& e = iter->second.front();
                _eweight[e] += w;
649 650
                _rec[e] += ecc;
                _drec[e] += decc;
651 652 653 654
            }
            else
            {
                auto e = add_edge(v, t, _g).first;
655
                ns_v[std::make_tuple(t, l)].push_back(e);
656
                eweight_c[e] = w;
657 658
                rec_c[e] = ecc;
                drec_c[e] = decc;
659 660 661 662 663 664 665 666 667 668
                set_prop(ec, e, l);
            }
        }

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

            for(auto e : in_edges_range(v, _g))
669
                ns_v[std::make_tuple(source(e, _g), ec[e])].push_back(e);
670
            for(auto e : in_edges_range(u, _g))
671
                ns_u[std::make_tuple(source(e, _g), ec[e])].push_back(e);
672 673 674 675 676 677 678 679 680 681 682

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

                if (s == u)
                    continue;

                size_t w = 0;
683
                std::vector<double> ecc, decc;
684
                for (auto& e : es)
685
                {
686
                    w += _eweight[e];
687 688 689
                    ecc += _rec[e];
                    decc += _drec[e];
                }
690

691
                auto iter = ns_v.find(std::make_tuple(s, l));
692 693 694 695
                if (iter != ns_v.end())
                {
                    auto& e = iter->second.front();
                    _eweight[e] += w;
696 697
                    _rec[e] += ecc;
                    _drec[e] += decc;
698 699 700 701
                }
                else
                {
                    auto e = add_edge(s, v, _g).first;
702
                    ns_v[std::make_tuple(s, l)].push_back(e);
703
                    eweight_c[e] = w;
704 705
                    rec_c[e] = ecc;
                    drec_c[e] = decc;
706 707 708 709 710 711 712 713
                    set_prop(ec, e, l);
                }
            }
        }

        _vweight[v] +=_vweight[u];
        _vweight[u] = 0;
        for (auto e : all_edges_range(u, _g))
714
        {
715
            _eweight[e] = 0;
716 717
            _rec[e].clear();
            _drec[e].clear();
718
        }
719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740
        clear_vertex(u, _g);
        _merge_map[u] = v;
        merge_degs(u, v, _degs);
    }

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

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

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

    void merge_degs(size_t u, size_t v, typename degs_map_t::unchecked_t& degs)
    {
        gt_hash_map<std::tuple<size_t, size_t>, size_t> hist;
        for (auto& kn : degs[u])
741
            hist[make_tuple(get<0>(kn), get<1>(kn))] += get<2>(kn);
742
        for (auto& kn : degs[v])
743
            hist[make_tuple(get<0>(kn), get<1>(kn))] += get<2>(kn);
744 745 746 747 748 749 750
        degs[u].clear();
        degs[v].clear();
        auto& d = degs[v];
        for (auto& kn : hist)
            d.emplace_back(get<0>(kn.first), get<1>(kn.first), kn.second);
    }

751 752 753
    // =========================================================================
    // Virtual state modification
    // =========================================================================
754 755 756 757 758 759 760

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

764
        double dS = entries_dS<exact>(m_entries, _mrs, _emat, _bg);
765

766
        size_t kout = out_degreeS()(v, _g, _eweight);
767 768
        size_t kin = kout;
        if (is_directed::apply<g_t>::type::value)
769
            kin = in_degreeS()(v, _g, _eweight);
770

771
        int dwr = _vweight[v];
772 773
        int dwnr = dwr;

774 775 776 777 778
        if (r == null_group && dwnr == 0)
            dwnr = 1;

        auto vt = [&](auto mrp, auto mrm, auto nr)
            {
779
                assert(mrp >= 0 && mrm >=0 && nr >= 0);
780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796
                if (exact)
                    return vterm_exact(mrp, mrm, nr, _deg_corr, _bg);
                else
                    return vterm(mrp, mrm, nr, _deg_corr, _bg);
            };

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

        if (nr != null_group)
        {
            dS += vt(_mrp[nr] + kout, _mrm[nr] + kin, _wr[nr] + dwnr);
            dS -= vt(_mrp[nr]       , _mrm[nr]      , _wr[nr]       );
        }
797 798 799 800

        return dS;
    }

801 802
    template <bool exact>
    double virtual_move_sparse(size_t v, size_t r, size_t nr)
803
    {
804
        return virtual_move_sparse<exact>(v, r, nr);
805 806
    }

807
    double virtual_move_dense(size_t v, size_t r, size_t nr, bool multigraph)
808 809 810 811 812 813 814 815 816 817
    {
        if (_deg_corr)
            throw GraphException("Dense entropy for degree corrected model not implemented!");

        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;

        if (r == nr)
            return 0;

        int kin = 0, kout = 0;
818
        kout += out_degreeS()(v, _g, _eweight);
819
        if (is_directed::apply<g_t>::type::value)
820
            kin += in_degreeS()(v, _g, _eweight);
821 822 823

        vector<int> deltap(num_vertices(_bg), 0);
        int deltal = 0;
824
        for (auto e : out_edges_range(v, _g))
825
        {
826 827
            vertex_t u = target(e, _g);
            vertex_t s = _b[u];
828
            if (u == v)
829
                deltal += _eweight[e];
830
            else
831
                deltap[s] += _eweight[e];
832 833 834 835 836
        }
        if (!is_directed::apply<g_t>::type::value)
            deltal /= 2;

        vector<int> deltam(num_vertices(_bg), 0);
837
        for (auto e : in_edges_range(v, _g))
838
        {
839
            vertex_t u = source(e, _g);
840 841
            if (u == v)
                continue;
842 843
            vertex_t s = _b[u];
            deltam[s] += _eweight[e];
844 845 846
        }

        double dS = 0;
847
        int dwr = _vweight[v];
848 849
        int dwnr = dwr;

850 851 852 853 854 855 856 857 858 859
        if (r == null_group && dwnr == 0)
            dwnr = 1;

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

860 861 862
        double Si = 0, Sf = 0;
        for (vertex_t s = 0; s < num_vertices(_bg); ++s)
        {
863 864
            int ers = (r != null_group) ? get_beprop(r, s, _mrs, _emat) : 0;
            int enrs = (nr != null_group) ? get_beprop(nr, s, _mrs, _emat) : 0;
865 866 867 868 869

            if (!is_directed::apply<g_t>::type::value)
            {
                if (s != nr && s != r)
                {
870 871 872 873 874 875 876 877 878 879 880
                    if (r != null_group)
                    {
                        Si += eterm_dense(r,  s, ers,              _wr[r],         _wr[s], multigraph, _bg);
                        Sf += eterm_dense(r,  s, ers - deltap[s],  _wr[r] - dwr,   _wr[s], multigraph, _bg);
                    }

                    if (nr != null_group)
                    {
                        Si += eterm_dense(nr, s, enrs,             _wr[nr],        _wr[s], multigraph, _bg);
                        Sf += eterm_dense(nr, s, enrs + deltap[s], _wr[nr] + dwnr, _wr[s], multigraph, _bg);
                    }
881 882 883 884 885 886 887 888 889 890 891 892 893
                }

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

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

894 895 896 897 898
                    if (r != null_group)
                    {
                        Si += eterm_dense(r, nr, ers,                          _wr[r],       _wr[nr],        multigraph, _bg);
                        Sf += eterm_dense(r, nr, ers - deltap[nr] + deltap[r], _wr[r] - dwr, _wr[nr] + dwnr, multigraph, _bg);
                    }
899 900 901 902
                }
            }
            else
            {
903 904
                int esr = (r != null_group) ? get_beprop(s, r, _mrs, _emat) : 0;
                int esnr  = (nr != null_group) ? get_beprop(s, nr, _mrs, _emat) : 0;
905 906 907

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

                    if (nr != null_group)
                    {
                        Si += eterm_dense(nr, s, enrs            , _wr[nr]       , _wr[s]        , multigraph, _bg);
                        Sf += eterm_dense(nr, s, enrs + deltap[s], _wr[nr] + dwnr, _wr[s]        , multigraph, _bg);
                        Si += eterm_dense(s, nr, esnr            , _wr[s]        , _wr[nr]       , multigraph, _bg);
                        Sf += eterm_dense(s, nr, esnr + deltam[s], _wr[s]        , _wr[nr] + dwnr, multigraph, _bg);
                    }
923 924 925 926 927 928 929
                }

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

930 931 932 933 934
                    if (nr != null_group)
                    {
                        Si += eterm_dense(r, nr, esnr                         , _wr[r]      , _wr[nr]       , multigraph, _bg);
                        Sf += eterm_dense(r, nr, esnr - deltap[nr] + deltam[r], _wr[r] - dwr, _wr[nr] + dwnr, multigraph, _bg);
                    }
935 936 937 938 939 940 941
                }

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

942 943 944 945 946
                    if (r != null_group)
                    {
                        Si += eterm_dense(nr, r, esr                         , _wr[nr]       , _wr[r]      , multigraph, _bg);
                        Sf += eterm_dense(nr, r, esr + deltap[r] - deltam[nr], _wr[nr] + dwnr, _wr[r] - dwr, multigraph, _bg);
                    }
947 948 949 950 951 952 953
                }
            }
        }

        return Sf - Si + dS;
    }

954
    template <class MEntries>
955
    double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea,
956
                        MEntries& m_entries)
957
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
958
        assert(size_t(_b[v]) == r || r == null_group);
959

960 961
        if (r == nr)
            return 0;
962