graph_blockmodel.hh

// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_BLOCKMODEL_HH
#define GRAPH_BLOCKMODEL_HH

#include "config.h"

#include <vector>

#include "../support/graph_state.hh"
#include "graph_blockmodel_util.hh"

#include "openmp_lock.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;

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

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;

#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

#define BLOCK_STATE_params                                                     \
    ((g, &, all_graph_views, 1))                                               \
    ((is_weighted,, bool_tr, 1))                                               \
    ((use_hash,, bool_tr, 1))                                                  \
    ((use_rmap,, rmap_tr, 1))                                                  \
    ((_abg, &, boost::any&, 0))                                                \
    ((_aeweight, &, boost::any&, 0))                                           \
    ((_avweight, &, boost::any&, 0))                                           \
    ((_adegs, &, boost::any&, 0))                                              \
    ((mrs,, emap_t, 0))                                                        \
    ((mrp,, vmap_t, 0))                                                        \
    ((mrm,, vmap_t, 0))                                                        \
    ((wr,, vmap_t, 0))                                                         \
    ((b,, vmap_t, 0))                                                          \
    ((empty_blocks, & ,std::vector<size_t>&, 0))                               \
    ((empty_pos,, vmap_t, 0))                                                  \
    ((candidate_blocks, &, std::vector<size_t>&, 0))                           \
    ((candidate_pos,, vmap_t, 0))                                              \
    ((bclabel,, vmap_t, 0))                                                    \
    ((pclabel,, vmap_t, 0))                                                    \
    ((bfield,, vprop_map_t<std::vector<double>>::type, 0))                     \
    ((Bfield, &, std::vector<double>&, 0))                                     \
    ((merge_map,, vmap_t, 0))                                                  \
    ((deg_corr,, bool, 0))                                                     \
    ((rec_types,, std::vector<int32_t>, 0))                                    \
    ((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))                      \
    ((brecsum,, vprop_map_t<double>::type, 0))                                 \
    ((wparams, &, std::vector<std::vector<double>>&, 0))                       \
    ((recdx, &, std::vector<double>&, 0))                                      \
    ((Lrecdx, &, std::vector<double>&, 0))                                     \
    ((epsilon, &, std::vector<double>&, 0))                                    \
    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))                \
    ((bignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))               \
    ((allow_empty,, bool, 0))

GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)

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

    typedef partition_stats<use_rmap_t::value> partition_stats_t;

    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()),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
          _emat(_bg, rng),
          _neighbor_sampler(_g, _eweight),
          _m_entries(num_vertices(_bg))
    {
        _empty_blocks.clear();
        _candidate_blocks.clear();
        _candidate_blocks.push_back(null_group);
        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);
        }
        for (auto& p : _rec)
            _c_rec.push_back(p.get_checked());
        for (auto& p : _drec)
            _c_drec.push_back(p.get_checked());
        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());

        if (!_rec_types.empty())
        {
            _recx2.resize(_rec_types.size());
            _recdx.resize(_rec_types.size());
            for (auto me : edges_range(_bg))
            {
                if (_brec[0][me] > 0)
                {
                    _B_E++;
                    for (size_t i = 0; i < _rec_types.size(); ++i)
                    {
                        if (_rec_types[i] == weight_type::REAL_NORMAL)
                        {
                            _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++;
            }
        }

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

    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()),
          _c_rec(other._c_rec),
          _c_drec(other._c_drec),
          _c_brec(other._c_brec),
          _c_bdrec(other._c_bdrec),
          _recsum(other._recsum),
          _recx2(other._recx2),
          _dBdx(other._dBdx),
          _LdBdx(other._LdBdx),
          _B_E(other._B_E),
          _B_E_D(other._B_E_D),
          _rt(other._rt),
          _vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
          _eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
          _degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
          _emat(other._emat),
          _egroups_enabled(other._egroups_enabled),
          _neighbor_sampler(other._neighbor_sampler),
          _m_entries(num_vertices(_bg))
    {
        if (other.is_partition_stats_enabled())
            enable_partition_stats();
    }

    // =========================================================================
    // State modification
    // =========================================================================

    template <class MEntries, class EFilt>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries,
                          EFilt&& efilt)
    {
        auto mv_entries = [&](auto&&... args)
            {
                move_entries(v, r, nr, _b, _g, _eweight, num_vertices(_bg),
                             m_entries, std::forward<EFilt>(efilt),
                             is_loop_nop(),
                             std::forward<decltype(args)>(args)...);
            };

        if (_rt == weight_type::NONE)
        {
            mv_entries();
        }
        else
        {
            if (_rt == weight_type::REAL_NORMAL)
                mv_entries(_rec, _drec);
            else
                mv_entries(_rec);
        }
    }

    template <class MEntries>
    void get_move_entries(size_t v, size_t r, size_t nr, MEntries& m_entries)
    {
        get_move_entries(v, r, nr, m_entries, [](auto) {return false;});
    }


    template <bool Add, class EFilt>
    void modify_vertex(size_t v, size_t r, EFilt&& efilt)
    {
        if (Add)
            get_move_entries(v, null_group, r, _m_entries,
                             std::forward<EFilt>(efilt));
        else
            get_move_entries(v, r, null_group, _m_entries,
                             std::forward<EFilt>(efilt));

        apply_delta<Add,!Add>(*this, _m_entries);

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

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

    bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true)
    {
        if (_coupled_state != nullptr && is_last(v))
        {
            auto& bh = _coupled_state->get_b();
            if (bh[r] != bh[nr])
                return false;
        }

        if (allow_empty)
            return ((_bclabel[r] == _bclabel[nr]) || (_wr[nr] == 0));
        else
            return _bclabel[r] == _bclabel[nr];
    }

    template <class EFilt>
    void move_vertex(size_t v, size_t r, size_t nr, EFilt&& efilt)
    {
        if (r == nr)
            return;

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

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

        apply_delta<true, true>(*this, _m_entries);

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

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

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

    void move_vertex(size_t v, size_t nr)
    {
        size_t r = _b[v];
        move_vertex(v, r, nr);
    }

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

        if (_rt == weight_type::NONE)
        {
            for (auto& rsd : entries)
            {
                _m_entries.template insert_delta<true>(_b[get<0>(rsd)],
                                                       _b[get<1>(rsd)],
                                                       get<3>(rsd));
                update_edge(get<2>(rsd));
            }
        }
        else
        {
            for (auto& rsd : entries)
            {
                recs_propagate_insert(*this, _b[get<0>(rsd)], _b[get<1>(rsd)],
                                      get<2>(rsd), get<3>(rsd), get<4>(rsd),
                                      _m_entries);
                update_edge(get<2>(rsd));
            }
        }
        apply_delta<true, true>(*this, _m_entries);
    }

    void add_edge(const GraphInterface::edge_t& e)
    {
        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)];
        auto me = _emat.get_me(r, s);
        if (me != _emat.get_null_edge() && _mrs[me] == 0)
        {
            _emat.remove_me(me, _bg);
            if (_coupled_state != nullptr)
                _coupled_state->remove_edge(me);
        }
        assert(e != _emat.get_null_edge());
        update_edge(e);
        boost::remove_edge(e, _g);
    }

    void update_edge(const GraphInterface::edge_t& e)
    {
        update_edge(e, is_weighted_t());
    }

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

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

    void add_edge_rec(const GraphInterface::edge_t& e)
    {
        if (_rec_types.empty())
            return;
        auto crec = _rec[0].get_checked();
        crec[e] = 1;
        for (size_t i = 1; i < _rec_types.size(); ++i)
        {
            auto drec = _drec[i].get_checked();
            drec[e] = 0;
        }
    }

    void remove_edge_rec(const GraphInterface::edge_t& e)
    {
        if (_rec_types.empty())
            return;
        _rec[0][e] = 0;
    }

    void update_edge_rec(const GraphInterface::edge_t& e,
                         const std::vector<double>& delta)
    {
        if (_rec_types.empty())
            return;

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

            auto rec = _c_rec[i][e];
            auto d = (std::pow(rec, 2) -
                      std::pow(rec - delta[i], 2));
            _c_drec[i][e] += d;
        }
    }

    void remove_partition_node(size_t v, size_t r)
    {
        assert(size_t(_b[v]) == r);

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

        _wr[r] -= _vweight[v];

        if (!_egroups.empty() && _egroups_enabled)
            _egroups.remove_vertex(v, _b, _g);

        if (is_partition_stats_enabled())
            get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g,
                                                 _vweight, _eweight,
                                                 _degs);
    }

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

        _wr[r] += _vweight[v];

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

        if (is_partition_stats_enabled())
            get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
                                              _eweight, _degs);

        if (_vweight[v] > 0 && _wr[r] == _vweight[v])
        {
            remove_element(_empty_blocks, _empty_pos, r);
            add_element(_candidate_blocks, _candidate_pos, r);

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

    template <class EFilt>
    void remove_vertex(size_t v, size_t r, EFilt&& efilt)
    {
        modify_vertex<false>(v, r, std::forward<EFilt>(efilt));
    }

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

    void remove_vertex(size_t v)
    {
        size_t r = _b[v];
        remove_vertex(v, r);
    }

    template <class Vlist>
    void remove_vertices(Vlist& vs)
    {
        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
        typedef typename graph_traits<g_t>::edge_descriptor edges_t;

        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)
            remove_vertex(v, _b[v],
                          [&](auto& e) { return eset.find(e) != eset.end(); });

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

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

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

            assert(_mrs[me] >= 0);

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

            for (size_t i = 0; i < _rec_types.size(); ++i)
            {
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
                    _bdrec[i][me] -= _drec[i][e];
                    [[gnu::fallthrough]];
                default:
                    _brec[i][me] -= _rec[i][e];
                }
            }

            if (_mrs[me] == 0)
            {
                _emat.remove_me(me, _bg);
                if (_coupled_state != nullptr)
                    _coupled_state->remove_edge(me);
                else
                    boost::remove_edge(me, this->_bg);
            }
        }
    }

    void remove_vertices(python::object ovs)
    {
        multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(ovs);
        remove_vertices(vs);
    }

    template <class EFilt>
    void add_vertex(size_t v, size_t r, EFilt&& efilt)
    {
        modify_vertex<true>(v, r, std::forward<EFilt>(efilt));
    }

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

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

        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)
            add_vertex(vr.first, vr.second,
                       [&](auto& e){ return eset.find(e) != eset.end(); });

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

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

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

            auto ew = _eweight[e];

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

            for (size_t i = 0; i < _rec_types.size(); ++i)
            {
                switch (_rec_types[i])
                {
                case weight_type::REAL_NORMAL: // signed weights
                    _bdrec[i][me] += _drec[i][e];
                    [[gnu::fallthrough]];
                default:
                    _brec[i][me] += _rec[i][e];
                }
            }
        }
    }

    void add_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);
        add_vertices(vs, rs);
    }

    template <bool Add, bool Deplete=true>
    void modify_edge(size_t u, size_t v, GraphInterface::edge_t& e,
                     const std::vector<double>& rec)
    {
        size_t r = _b[u];
        size_t s = _b[v];

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

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

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

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

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

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

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

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

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

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

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

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

    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));
        vweight[v] = 0;
    }

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

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

    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)
    {
        merge_vertices(u, v, ec, _is_weighted);
    }

    template <class Emap>
    void merge_vertices(size_t, size_t, Emap&&, std::false_type)
    {
        throw ValueException("cannot merge vertices of unweighted graph");
    }

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

        auto eweight_c = _eweight.get_checked();
        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());

        typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
        typedef typename graph_traits<g_t>::edge_descriptor edge_t;

        gt_hash_map<std::tuple<vertex_t, int>, vector<edge_t>> ns_u, ns_v;
        for(auto e : out_edges_range(u, _g))
            ns_u[std::make_tuple(target(e, _g), ec[e])].push_back(e);
        for(auto e : out_edges_range(v, _g))
            ns_v[std::make_tuple(target(e, _g), ec[e])].push_back(e);

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

        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;

            std::fill(ecc.begin(), ecc.end(), 0);
            std::fill(decc.begin(), decc.end(), 0);
            for (auto& e : es)
            {