graph_blockmodel.hh

// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2015 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 <cmath>
#include <iostream>
#include <queue>
#include <boost/math/special_functions/zeta.hpp>
#include <boost/functional/hash.hpp>

#include "config.h"
#include <unordered_set>
#include <unordered_map>
#include <tuple>

#ifdef HAVE_SPARSEHASH
#include SPARSEHASH_INCLUDE(dense_hash_set)
#include SPARSEHASH_INCLUDE(dense_hash_map)
#endif

#include "../generation/sampler.hh"
#include "../generation/dynamic_sampler.hh"

#ifdef USING_OPENMP
#include <omp.h>
#endif


double spence(double);

namespace graph_tool
{

#ifdef HAVE_SPARSEHASH
using google::dense_hash_set;
using google::dense_hash_map;
#endif

using namespace boost;

template <class Key>
class IdentityArrayPropertyMap
    : public boost::put_get_helper<Key, IdentityArrayPropertyMap<Key>>
{
public:
    typedef std::array<Key, 1> value_type;
    typedef value_type reference;
    typedef Key key_type;
    typedef boost::readable_property_map_tag category;

    inline __attribute__((always_inline))
    const value_type operator[](const key_type& c) const { return {c}; }
};

// ====================
// Entropy calculation
// ====================

// Repeated computation of x*log(x) and log(x) actually adds up to a lot of
// time. A significant speedup can be made by caching pre-computed values. This
// is doable since the values of mrse are bounded in [0, 2E], where E is the
// total number of edges in the network. Here we simply grow the cache as
// needed.

extern vector<double> __safelog_cache;
extern vector<double> __xlogx_cache;
extern vector<double> __lgamma_cache;

template <class Type>
__attribute__((always_inline))
inline double safelog(Type x)
{
    if (x == 0)
        return 0;
    return log(x);
}

__attribute__((always_inline))
inline double safelog(size_t x)
{
    assert(x < __safelog_cache.size());
    return __safelog_cache[x];
}

__attribute__((always_inline))
inline double xlogx(size_t x)
{
    if (x >= __xlogx_cache.size())
        cout << x << " " << __xlogx_cache.size() << endl;
    assert(x < __xlogx_cache.size());
    return __xlogx_cache[x];
}

__attribute__((always_inline))
inline double lgamma_fast(size_t x)
{
    if (x >= __lgamma_cache.size())
        return lgamma(x);
    return __lgamma_cache[x];
}

// polylogarithm and degree-distribution description length (xi)

template <class Type>
Type polylog(int n, Type z, Type epsilon=1e-6)
{
    if (n == 2)
        return spence(1 - z);

    int k = 1;
    Type S = 0;
    Type delta = epsilon + 1;
    Type zk = z;
    while (delta > epsilon)
    {
        Type dS = zk / pow(k, n);
        k++;
        zk *= z;
        S += dS;
        delta = dS;
    }
    return S;
}

template <class NType, class EType>
void get_mu_l(NType N, EType E, double& mu, double& l,
              double epsilon=1e-8)
{
    mu = sqrt(polylog<double>(2, 1.) / double(E));
    l = 1. - exp(-double(N) * mu);

    double delta = epsilon + 1;
    while (delta > epsilon)
    {
        double nmu = sqrt(polylog<double>(2, l) / double(E));
        double nl = 1. - exp(-double(N) * mu);

        delta = abs(nmu - mu) + abs(nl - l);
        mu = nmu;
        l = nl;
    }

    l = -log(l);
}

template <class NType, class EType>
double get_xi(NType N, EType E, double epsilon=1e-8)
{
    if (E == 0 || N == 0)
        return 0;

    double mu = 0, l = 0;
    get_mu_l(N, E, mu, l, epsilon);
    double S = double(N) * l + 2 * double(E) * mu;
    return S;
}

template <class NType, class EType>
double get_xi_fast(NType N, EType E)
{
    if (E == 0 || N == 0)
        return 0;
    static const double z2 = boost::math::zeta(2.);
    return 2 * sqrt(z2 * E);
}

// Sparse entropy terms
// ====================

// "edge" term of the entropy
template <class Graph>
__attribute__((always_inline))
inline double eterm(size_t r, size_t s, size_t mrs, const Graph&)
{
    if (!is_directed::apply<Graph>::type::value && r == s)
        mrs *= 2;

    double val = xlogx(mrs);

    if (is_directed::apply<Graph>::type::value || r != s)
        return -val;
    else
        return -val / 2;
}

// "vertex" term of the entropy
template <class Graph>
inline double vterm(size_t mrp, size_t mrm, size_t wr, bool deg_corr,
                    Graph&)
{
    double one = 0.5;

    if (is_directed::apply<Graph>::type::value)
        one = 1;

    if (deg_corr)
        return one * (xlogx(mrm) + xlogx(mrp));
    else
        return one * (mrm * safelog(wr) + mrp * safelog(wr));
}

struct entropy
{
    template <class Graph, class Eprop, class Vprop>
    void operator()(Eprop& mrs, Vprop& mrp, Vprop& mrm, Vprop& wr,
                    bool deg_corr, Graph& g, double& S) const
    {
        S = 0;
        for (auto e : edges_range(g))
            S += eterm(source(e, g), target(e, g), mrs[e], g);
        for (auto v : vertices_range(g))
            S += vterm(mrp[v], mrm[v], wr[v], deg_corr, g);
    }
};

// Parallel edges
// ==============

template <class Vertex, class List, class Graph, class GetNode>
double get_parallel_neighbours_entropy(Vertex v, List& us, Graph&,
                                       const GetNode& get_node)
{
    double S = 0;
    for (auto& uc : us)
    {
        auto& u = uc.first;
        auto& m = uc.second;
        if (m > 1)
        {
            if (get_node(u) == size_t(v) && !is_directed::apply<Graph>::type::value)
            {
                assert(m % 2 == 0);
                S += lgamma_fast(m/2 + 1);
            }
            else
            {
                S += lgamma_fast(m + 1);
            }
        }
    }
    return S;
}

struct entropy_parallel_edges
{
    template <class Graph, class Weight>
    void operator()(Graph& g, Weight weight, double& S) const
    {
        S = 0;
        auto get_node = [](size_t i) {return i;};
        for (auto v : vertices_range(g))
        {
            unordered_map<decltype(v), int> us;
            for (auto e : out_edges_range(v, g))
            {
                auto u = target(e, g);
                if (u < v && !is_directed::apply<Graph>::type::value)
                    continue;
                us[u] += weight[e];
            }

            S += get_parallel_neighbours_entropy(v, us, g, get_node);
        }
    }
};


// Dense entropy
// =============

__attribute__((always_inline))
inline double lbinom(double N, double k)
{
    if (N == 0 || k == 0 || k > N)
        return 0;
    return lgamma(N + 1) - lgamma(N - k + 1) - lgamma(k + 1);
}

__attribute__((always_inline))
inline double lbinom_fast(int N, int k)
{
    if (N == 0 || k == 0 || k > N)
        return 0;
    return lgamma_fast(N + 1) - lgamma_fast(N - k + 1) - lgamma_fast(k + 1);
}


// "edge" term of the entropy
template <class Graph>
__attribute__((always_inline))
inline double eterm_dense(size_t r, size_t s, int ers, double wr_r,
                          double wr_s, bool multigraph, const Graph&)
{
    // we should not use integers here, since they may overflow
    double nrns;

    if (ers == 0)
        return 0.;

    if (r != s || is_directed::apply<Graph>::type::value)
    {
        nrns = wr_r * wr_s;
    }
    else
    {
        if (multigraph)
            nrns = (wr_r * (wr_r + 1)) / 2;
        else
            nrns = (wr_r * (wr_r - 1)) / 2;
    }

    double S;
    if (multigraph)
        S = lbinom(nrns + ers - 1, ers); // do not use lbinom_fast!
    else
        S = lbinom(nrns, ers);
    return S;
}

struct entropy_dense
{
    template <class Graph, class Eprop, class Vprop>
    void operator()(Eprop mrs, Vprop& wr, bool multigraph, Graph& g, double& S) const
    {
        S = 0;
        for (auto e : edges_range(g))
        {
            auto r = source(e, g);
            auto s = target(e, g);
            S += eterm_dense(r, s, mrs[e], wr[r], wr[s], multigraph, g);
        }
    }
};

// ===============
// Partition stats
// ===============

class partition_stats_t
{
public:

#ifdef HAVE_SPARSEHASH
    typedef dense_hash_map<pair<size_t,size_t>, int, std::hash<pair<size_t,size_t>>> map_t;
#else
    typedef unordered_map<pair<size_t,size_t>, int> map_t;
#endif

    partition_stats_t() : _enabled(false) {}

    template <class Graph, class Vprop, class Eprop>
    partition_stats_t(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
                      bool edges_dl)
        : _enabled(true), _N(N), _E(0), _B(B), _hist(B), _total(B), _ep(B),
          _em(B), _edges_dl(edges_dl)
    {

#ifdef HAVE_SPARSEHASH
        for (size_t r = 0; r < B; ++r)
        {
            _hist[r].set_empty_key(make_pair(numeric_limits<size_t>::max(),
                                             numeric_limits<size_t>::max()));
            _hist[r].set_deleted_key(make_pair(numeric_limits<size_t>::max() - 1,
                                               numeric_limits<size_t>::max() - 1));
        }
#endif

        for (auto v : vertices_range(g))
        {
            auto r = b[v];
            size_t kin = in_degreeS()(v, g, eweight);
            size_t kout = out_degreeS()(v, g, eweight);
            _hist[r][make_pair(kin, kout)]++;
            _total[r]++;
            _em[r] += kin;
            _ep[r] += kout;
            _E += kout;
        }

        if (!is_directed::apply<Graph>::type::value)
            _E /= 2;

        _actual_B = 0;
        for (size_t r = 0; r < B; ++r)
            if (_total[r] > 0)
                _actual_B++;
    }

    double get_partition_dl()
    {
        double S = 0;
        S += lbinom(_actual_B + _N - 1, _N);
        S += lgamma(_N + 1);
        for (auto nr : _total)
            S -= lgamma(nr + 1);
        return S;
    }

    double get_deg_dl(bool ent, bool dl_alt, bool xi_fast)
    {
        double S = 0;
        for (size_t r = 0; r < _B; ++r)
        {
            if (ent)
            {
                for (auto& k_c : _hist[r])
                {
                    double p = k_c.second / double(_total[r]);
                    S -= p * log(p) * _total[r];
                }
            }
            else
            {
                double S1 = 0;

                if (xi_fast)
                {
                    S1 += get_xi_fast(_total[r], _ep[r]);
                    S1 += get_xi_fast(_total[r], _em[r]);
                }
                else
                {
                    S1 += get_xi(_total[r], _ep[r]);
                    S1 += get_xi(_total[r], _em[r]);
                }

                S1 += lgamma(_total[r] + 1);
                for (auto& k_c : _hist[r])
                    S1 -= lgamma(k_c.second + 1);

                if (dl_alt)
                {
                    double S2 = 0;
                    S2 += lbinom(_total[r] + _ep[r] - 1, _ep[r]);
                    S2 += lbinom(_total[r] + _em[r] - 1, _em[r]);
                    S += min(S1, S2);
                }
                else
                {
                    S += S1;
                }
            }
        }
        return S;
    }

    template <class Graph, class OStats>
    double get_delta_dl(size_t, size_t r, size_t nr, OStats&, Graph&)
    {
        if (r == nr)
            return 0;

        double S_b = 0, S_a = 0;
        S_b += -lgamma_fast(_total[r] + 1) - lgamma_fast(_total[nr] + 1);
        S_a += -lgamma_fast(_total[r]    ) - lgamma_fast(_total[nr] + 2);

        int dB = 0;
        if (_total[r] == 1)
            dB--;
        if (_total[nr] == 0)
            dB++;

        if (dB != 0)
        {
            S_b += lbinom(_actual_B + _N - 1, _N);
            S_a += lbinom(_actual_B + dB + _N - 1, _N);

            if (_edges_dl)
            {
                auto get_x = [](size_t B) -> size_t
                {
                    if (is_directed::apply<Graph>::type::value)
                        return B * B;
                    else
                        return (B * (B + 1)) / 2;
                };

                S_b += lbinom(get_x(_actual_B) + _E - 1, _E);
                S_a += lbinom(get_x(_actual_B + dB) + _E - 1, _E);
            }
        }

        return S_a - S_b;
    }

    template <class Graph, class EWeight, class OStats>
    double get_delta_deg_dl(size_t v, size_t r, size_t nr, EWeight& eweight,
                            OStats&, Graph& g)
    {
        if (r == nr)
            return 0;

        double S_b = 0, S_a = 0;

        int kin = in_degreeS()(v, g, eweight);
        int kout = out_degreeS()(v, g, eweight);

        auto get_Se = [&](size_t s, int delta, int kin, int kout) -> double
            {
                double S = 0;
                S += get_xi_fast(_total[s] + delta, _em[s] + kin);
                S += get_xi_fast(_total[s] + delta, _ep[s] + kout);
                return S;
            };

        S_b += get_Se(r,  0,    0,     0) + get_Se(nr, 0,   0,    0);
        S_a += get_Se(r, -1, -kin, -kout) + get_Se(nr, 1, kin, kout);

        auto get_Sr = [&](size_t s, int delta) -> double
            {
                return lgamma_fast(_total[s] + delta + 1);
            };

        S_b += get_Sr(r,  0) + get_Sr(nr, 0);
        S_a += get_Sr(r, -1) + get_Sr(nr, 1);


        auto get_Sk = [&](size_t s, pair<size_t, size_t>& deg, int delta) -> double
            {
                size_t nd = 0;
                auto iter = _hist[s].find(deg);
                if (iter != _hist[s].end())
                    nd = iter->second;

                return -lgamma_fast(nd + delta + 1);
            };

        auto deg = make_pair(size_t(kin), size_t(kout));
        S_b += get_Sk(r, deg,  0) + get_Sk(nr, deg, 0);
        S_a += get_Sk(r, deg, -1) + get_Sk(nr, deg, 1);

        return S_a - S_b;
    }

    template <class Graph, class OStats>
    void move_vertex(size_t v, size_t r, size_t nr, bool deg_corr, OStats&,
                     Graph& g, size_t kin = 0, size_t kout = 0)
    {
        if (r == nr)
            return;

        _total[r]--;
        _total[nr]++;

        if (_total[r] == 0)
            _actual_B--;
        if (_total[nr] == 1)
            _actual_B++;

        if (deg_corr)
        {
            if (kin + kout == 0)
            {
                kin = in_degreeS()(v, g);
                kout = out_degreeS()(v, g);
            }
            auto deg = make_pair(kin, kout);
            auto iter = _hist[r].find(deg);
            iter->second--;
            if (iter->second == 0)
                _hist[r].erase(iter);
            _hist[nr][deg]++;
            _em[r] -= deg.first;
            _ep[r] -= deg.second;
            _em[nr] += deg.first;
            _ep[nr] += deg.second;
        }
    }

    bool is_enabled() { return _enabled; }
    void set_enabled(bool enabled) { _enabled = enabled; }

private:
    bool _enabled;
    size_t _N;
    size_t _E;
    size_t _B;
    size_t _actual_B;
    vector<map_t> _hist;
    vector<int> _total;
    vector<int> _ep;
    vector<int> _em;
    bool _edges_dl;
};

// ===============================
// Block moves
// ===============================

// this structure speeds up the access to the edges between given blocks,
// since we're using an adjacency list to store the block structure (the emat_t
// is simply a corresponding adjacency matrix)
struct get_emat_t
{
    template <class Graph>
    struct apply
    {
        typedef multi_array<pair<typename graph_traits<Graph>::edge_descriptor, bool>, 2> type;
    };
};


struct create_emat
{
    template <class Graph>
    void operator()(Graph& g, boost::any& oemap) const
    {
        typedef typename get_emat_t::apply<Graph>::type emat_t;
        size_t B = num_vertices(g);
        emat_t emat(boost::extents[B][B]);

        for (size_t i = 0; i < B; ++i)
            for (size_t j = 0; j < B; ++j)
                emat[i][j].second = false;
        for (auto e : edges_range(g))
        {
            if (source(e, g) >= B || target(e, g) >= B)
                throw GraphException("incorrect number of blocks when creating emat!");
            emat[source(e, g)][target(e, g)] = make_pair(e, true);
            if (!is_directed::apply<Graph>::type::value)
                emat[target(e, g)][source(e, g)] = make_pair(e, true);
        }

        oemap = emat;
    }
};

template <class Graph>
inline __attribute__((always_inline))
pair<typename graph_traits<Graph>::edge_descriptor, bool>
get_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename get_emat_t::apply<Graph>::type& emat, const Graph&)
{
    return emat[r][s];
}

template <class Graph>
inline __attribute__((always_inline))
void
put_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename graph_traits<Graph>::edge_descriptor& e,
       typename get_emat_t::apply<Graph>::type& emat, const Graph&)
{
    emat[r][s] = make_pair(e, true);
    if (!is_directed::apply<Graph>::type::value && r != s)
        emat[s][r] = make_pair(e, true);
}

template <class Graph>
inline __attribute__((always_inline))
void
remove_me(typename graph_traits<Graph>::vertex_descriptor r,
          typename graph_traits<Graph>::vertex_descriptor s,
          const typename graph_traits<Graph>::edge_descriptor&,
          typename get_emat_t::apply<Graph>::type& emat, Graph&,
          bool delete_edge=true)
{
    if (!delete_edge)
    {
        emat[r][s].second = false;
        if (!is_directed::apply<Graph>::type::value && r != s)
            emat[s][r].second = false;
    }
}


// this structure speeds up the access to the edges between given blocks, since
// we're using an adjacency list to store the block structure (this is like
// emat_t above, but takes less space and is slower)
struct get_ehash_t
{
    template <class Graph>
    struct apply
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
#ifdef HAVE_SPARSEHASH
        typedef dense_hash_map<vertex_t, edge_t, std::hash<vertex_t>> map_t;
#else
        typedef unordered_map<vertex_t, edge_t> map_t;
#endif
        typedef vector<map_t> type;
    };
};


template<class Graph>
inline __attribute__((always_inline))
pair<typename graph_traits<Graph>::edge_descriptor, bool>
get_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename get_ehash_t::apply<Graph>::type& ehash, const Graph&)
{
    assert(r < ehash.size());
    const auto& map = ehash[r];
    auto iter = map.find(s);
    if (iter == map.end())
        return (make_pair(typename graph_traits<Graph>::edge_descriptor(), false));
    return make_pair(iter->second, true);
}

template<class Graph>
inline __attribute__((always_inline))
void
put_me(typename graph_traits<Graph>::vertex_descriptor r,
       typename graph_traits<Graph>::vertex_descriptor s,
       const typename graph_traits<Graph>::edge_descriptor& e,
       typename get_ehash_t::apply<Graph>::type& ehash,
       const Graph&)
{
    assert(r < ehash.size());
    ehash[r][s] = e;
    if (!is_directed::apply<Graph>::type::value)
        ehash[s][r] = e;
}

template<class Graph>
inline __attribute__((always_inline))
void
remove_me(typename graph_traits<Graph>::vertex_descriptor r,
          typename graph_traits<Graph>::vertex_descriptor s,
          const typename graph_traits<Graph>::edge_descriptor& e,
          typename get_ehash_t::apply<Graph>::type& ehash, Graph& bg,
          bool delete_edge=true)
{
    assert(r < ehash.size());
    ehash[r].erase(s);
    if (!is_directed::apply<Graph>::type::value)
        ehash[s].erase(r);
    if (delete_edge)
        remove_edge(e, bg);
}

struct create_ehash
{
    template <class Graph>
    void operator()(Graph& g, boost::any& oemap) const
    {
        typedef typename get_ehash_t::apply<Graph>::type emat_t;
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

        emat_t emat(num_vertices(g));

#ifdef HAVE_SPARSEHASH
        for (auto v : vertices_range(g))
        {
            emat[v].set_empty_key(numeric_limits<vertex_t>::max());
            emat[v].set_deleted_key(numeric_limits<vertex_t>::max() - 1);
        }
#endif

        for (auto e : edges_range(g))
            put_me(source(e, g), target(e, g), e, emat, g);

#ifdef HAVE_SPARSEHASH
        for (auto v : vertices_range(g))
            emat[v].resize(0);
#endif
        oemap = emat;
    }
};

template <class Vertex, class Eprop, class Emat, class BGraph>
__attribute__((always_inline))
inline size_t get_mrs(Vertex r, Vertex s, const Eprop& mrs, Emat& emat,
                      BGraph& bg)
{
    const pair<typename graph_traits<BGraph>::edge_descriptor, bool> me =
        get_me(r, s, emat, bg);
    if (me.second)
        return mrs[me.first];
    else
        return 0;
}

struct standard_neighbours_policy
{
    template <class Graph, class Vertex>
    IterRange<typename out_edge_iteratorS<Graph>::type>
    get_out_edges(Vertex v, Graph& g) const
    {
        return out_edges_range(v, g);
    }

    template <class Graph, class Vertex>
    IterRange<typename in_edge_iteratorS<Graph>::type>
    get_in_edges(Vertex v, Graph& g) const
    {
        return in_edges_range(v, g);
    }

    template <class Graph, class Vertex, class Weight>
    int get_out_degree(Vertex& v, Graph& g, Weight& eweight) const
    {
        return out_degreeS()(v, g, eweight);
    }

    template <class Graph, class Vertex, class Weight>
    int get_in_degree(Vertex& v, Graph& g, Weight& eweight) const
    {
        return in_degreeS()(v, g, eweight);
    }
};

// remove a vertex from its current block
template <class Graph, class BGraph, class Eprop, class Vprop, class EWprop,
          class VWprop, class EMat, class OStats,
          class NPolicy = standard_neighbours_policy>
void remove_vertex(size_t v, Eprop& mrs, Vprop& mrp, Vprop& mrm, Vprop& wr,
                   Vprop& b, const EWprop& eweight, const VWprop& vweight,
                   Graph& g, BGraph& bg, EMat& emat, OStats& overlap_stats,
                   const NPolicy& npolicy = NPolicy())
{
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    vertex_t r = b[v];

    int self_weight = 0;
    for (auto e : npolicy.get_out_edges(v, g))
    {
        vertex_t u = target(e, g);
        vertex_t s = b[u];

        // if (!emat[r][s].second)
        //     throw GraphException("no edge? " + lexical_cast<string>(r) +
        //                          " " + lexical_cast<string>(s));

        auto me = get_me(r, s, emat, bg).first;

        size_t ew = eweight[e];
        if (u == v && !is_directed::apply<Graph>::type::value)
        {
            self_weight += ew;
        }
        else
        {
            mrs[me] -= ew;

            assert(mrs[me] >= 0);

            mrp[r] -= ew;
            mrm[s] -= ew;

            if (mrs[me] == 0)
                remove_me(r, s, me, emat, bg);
        }
    }

    if (self_weight > 0)
    {
        assert(self_weight % 2 == 0);
        auto me = get_me(r, r, emat, bg).first;
        mrs[me] -= self_weight / 2;
        mrp[r] -= self_weight / 2;
        mrm[r] -= self_weight / 2;
        assert(mrs[me] >= 0);
        if (mrs[me] == 0)
            remove_me(r, r, me, emat, bg);
    }

    for (auto e : npolicy.get_in_edges(v, g))
    {
        vertex_t u = source(e, g);
        if (u == v)
            continue;
        vertex_t s = b[u];

        // if (!emat[s][r].second)
        //     throw GraphException("no edge? " + lexical_cast<string>(s) +
        //                          " " + lexical_cast<string>(r));

        typename graph_traits<BGraph>::edge_descriptor me =
            get_me(s, r, emat, bg).first;

        size_t ew = eweight[e];
        mrs[me] -= ew;

        mrp[s] -= ew;
        mrm[r] -= ew;

        if (mrs[me] == 0)
            remove_me(s, r, me, emat, bg);
    }

    if (!overlap_stats.is_enabled())
    {
        wr[r] -= vweight[v];
    }
    else
    {
        overlap_stats.remove_half_edge(v, r, b, g);
        wr[r] = overlap_stats.get_block_size(r);
    }
}

// add a vertex to block rr
template <class Graph, class BGraph, class Eprop, class Vprop, class EWprop,
          class VWprop, class EMat, class OStats,
          class NPolicy = standard_neighbours_policy>
void add_vertex(size_t v, size_t r, Eprop& mrs, Vprop& mrp, Vprop& mrm,
                Vprop& wr, Vprop& b, const EWprop& eweight,
                const VWprop& vweight, Graph& g, BGraph& bg, EMat& emat,
                OStats& overlap_stats, const NPolicy& npolicy = NPolicy())
{
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;

    int self_weight = 0;
    for (auto e : npolicy.get_out_edges(v, g))
    {
        vertex_t u = target(e, g);
        vertex_t s;

        if (u != v)
            s = b[u];
        else
            s = r;

        typename graph_traits<BGraph>::edge_descriptor me;

        auto mep = get_me(r, s, emat, bg);

        if (!mep.second)
        {
            mep = add_edge(r, s, bg);
            put_me(r, s, mep.first, emat, bg);
            mrs[mep.first] = 0;
        }
        me = mep.first;

        size_t ew = eweight[e];

        if (u == v && !is_directed::apply<Graph>::type::value)
        {
            self_weight += ew;
        }
        else
        {
            mrs[me] += ew;
            mrp[r] += ew;
            mrm[s] += ew;
        }
    }

    if (self_weight > 0)
    {
        assert(self_weight % 2 == 0);
        auto me = get_me(r, r, emat, bg).first;
        mrs[me] += self_weight / 2;
        mrp[r] += self_weight / 2;
        mrm[r] += self_weight / 2;
        assert(mrs[me] >= 0);
    }

    for (auto e : npolicy.get_in_edges(v, g))
    {
        vertex_t u = source(e, g);
        if (u == v)
            continue;

        vertex_t s = b[u];

        typename graph_traits<BGraph>::edge_descriptor me;
        pair<typename graph_traits<BGraph>::edge_descriptor, bool> mep =
                get_me(s, r, emat, bg);


        if (!mep.second)
        {
            mep = add_edge(s, r, bg);
            put_me(s, r, mep.first, emat, bg);
            mrs[mep.first] = 0;
        }
        me = mep.first;

        size_t ew = eweight[e];

        mrs[me] += ew;

        mrp[s] += ew;
        mrm[r] += ew;
    }

    if (!overlap_stats.is_enabled())