Skip to content

GitLab

graph-tool
graph-tool
Commit 36514eb1 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files
Options

Fix bug in BlockState.entropy() when deg_corr=True and edge weights are being used 

This also adds the 'ignore_degrees' parameter.
parent 254b7f3a
Hide whitespace changes
Inline Side-by-side
Showing with 79 additions and 34 deletions
src/graph/community/graph_blockmodel.cc
View file @ 36514eb1
......@@ -585,13 +585,15 @@ void do_collect_vertex_marginals(GraphInterface& gi, boost::any ob,
struct get_deg_entropy_term
{
template <class Graph>
void operator()(Graph& g, double& S) const
template <class Graph, class Emap, class VImap>
void operator()(Graph& g, Emap eweight, VImap ignore_degrees, double& S) const
{
for(auto v : vertices_range(g))
{
S -= lgamma_fast(out_degree(v, g) + 1);
S -= lgamma_fast(in_degreeS()(v, g) + 1);
if (ignore_degrees[v])
continue;
S -= lgamma_fast(out_degreeS()(v, g, eweight) + 1);
S -= lgamma_fast(in_degreeS()(v, g, eweight) + 1);
}
}
};
......@@ -633,11 +635,18 @@ struct get_deg_entropy_term_overlap
};
double do_get_deg_entropy_term(GraphInterface& gi, boost::any ob,
overlap_stats_t& overlap_stats, size_t N)
overlap_stats_t& overlap_stats, size_t N,
boost::any aeweight, boost::any aignore_degrees)
{
typedef property_map_type::apply<int32_t,
GraphInterface::vertex_index_map_t>::type
vmap_t;
typedef property_map_type::apply<int32_t,
GraphInterface::edge_index_map_t>::type
emap_t;
typedef property_map_type::apply<uint8_t,
GraphInterface::vertex_index_map_t>::type
vimap_t;
double S = 0;
......@@ -651,9 +660,12 @@ double do_get_deg_entropy_term(GraphInterface& gi, boost::any ob,
}
else
{
emap_t eweight = any_cast<emap_t>(aeweight);
vimap_t ignore_degrees = any_cast<vimap_t>(aignore_degrees);
run_action<>()
(gi, std::bind(get_deg_entropy_term(),
placeholders::_1, std::ref(S)))();
placeholders::_1, eweight, ignore_degrees,
std::ref(S)))();
}
return S;
......@@ -721,17 +733,20 @@ python::tuple python_get_mu_l(double N, double E, double epsilon=1e-8)
struct get_partition_stats
{
template <class Graph, class Vprop, class Eprop>
template <class Graph, class Vprop, class Eprop, class Mprop>
void operator()(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
bool edges_dl, partition_stats_t& partition_stats) const
bool edges_dl, partition_stats_t& partition_stats,
Mprop ignore_degrees) const
{
partition_stats = partition_stats_t(g, b, eweight, N, B, edges_dl);
partition_stats = partition_stats_t(g, b, eweight, N, B, edges_dl,
ignore_degrees);
}
};
partition_stats_t
do_get_partition_stats(GraphInterface& gi, boost::any ob, boost::any aeweight,
size_t N, size_t B, bool edges_dl)
size_t N, size_t B, bool edges_dl,
boost::any aignore_degrees)
{
typedef property_map_type::apply<int32_t,
GraphInterface::vertex_index_map_t>::type
......@@ -739,15 +754,19 @@ do_get_partition_stats(GraphInterface& gi, boost::any ob, boost::any aeweight,
typedef property_map_type::apply<int32_t,
GraphInterface::edge_index_map_t>::type
emap_t;
typedef property_map_type::apply<uint8_t,
GraphInterface::vertex_index_map_t>::type
mvmap_t;
partition_stats_t partition_stats;
vmap_t b = any_cast<vmap_t>(ob);
emap_t eweight = any_cast<emap_t>(aeweight);
mvmap_t ignore_degrees = any_cast<mvmap_t>(aignore_degrees);
run_action<>()(gi, std::bind(get_partition_stats(),
placeholders::_1, b, eweight, N, B, edges_dl,
std::ref(partition_stats)))();
std::ref(partition_stats), ignore_degrees))();
return partition_stats;
}
......
src/graph/community/graph_blockmodel.hh
View file @ 36514eb1
......@@ -376,21 +376,27 @@ public:
partition_stats_t() : _enabled(false) {}
template <class Graph, class Vprop, class Eprop>
template <class Graph, class Vprop, class Eprop, class Mprop>
partition_stats_t(Graph& g, Vprop b, Eprop eweight, size_t N, size_t B,
bool edges_dl)
bool edges_dl, Mprop ignore_degree)
: _enabled(true), _N(N), _E(0), _B(B), _hist(B), _total(B), _ep(B),
_em(B), _edges_dl(edges_dl)
_em(B), _edges_dl(edges_dl), _ignore_degree(N, false)
{
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)]++;
if (v >= _ignore_degree.size())
_ignore_degree.resize(v + 1, false);
_ignore_degree[v] = ignore_degree[v];
if (!_ignore_degree[v])
{
_hist[r][make_pair(kin, kout)]++;
_em[r] += kin;
_ep[r] += kout;
}
_total[r]++;
_em[r] += kin;
_ep[r] += kout;
_E += kout;
}
......@@ -418,6 +424,9 @@ public:
double S = 0;
for (size_t r = 0; r < _B; ++r)
{
if (_ep[r] + _em[r] == 0)
continue;
if (ent)
{
for (auto& k_c : _hist[r])
......@@ -504,7 +513,7 @@ public:
double get_delta_deg_dl(size_t v, size_t r, size_t nr, EWeight& eweight,
OStats&, Graph& g)
{
if (r == nr)
if (r == nr || _ignore_degree[v])
return 0;
double S_b = 0, S_a = 0;
......@@ -531,7 +540,6 @@ public:
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;
......@@ -565,7 +573,7 @@ public:
if (_total[nr] == 1)
_actual_B++;
if (deg_corr)
if (deg_corr && !_ignore_degree[v])
{
if (kin + kout == 0)
{
......@@ -599,6 +607,7 @@ private:
vector<int> _ep;
vector<int> _em;
bool _edges_dl;
vector<bool> _ignore_degree;
};
// ===============================
......
src/graph_tool/community/blockmodel.py
View file @ 36514eb1
......@@ -194,6 +194,9 @@ class BlockState(object):
self.max_BE = max_BE
self.overlap = False
self.ignore_degrees = kwargs.get("ignore_degrees", None)
if self.ignore_degrees is None:
self.ignore_degrees = g.new_vertex_property("bool", False)
# used by mcmc_sweep()
self.egroups = None
......@@ -233,7 +236,8 @@ class BlockState(object):
_prop("v", self.g, self.b),
_prop("e", self.g, self.eweight),
self.N, self.B,
edges_dl)
edges_dl,
_prop("v", self.g, self.ignore_degrees))
else:
self.partition_stats = libcommunity.partition_stats()
......@@ -253,7 +257,8 @@ class BlockState(object):
B=(self.B if b is None else None) if B is None else B,
clabel=self.clabel if clabel is None else clabel,
deg_corr=self.deg_corr if deg_corr is None else deg_corr,
max_BE=self.max_BE)
max_BE=self.max_BE,
ignore_degrees=self.ignore_degrees)
else:
state = OverlapBlockState(self.g,
b=b if b is not None else self.b,
......@@ -281,7 +286,8 @@ class BlockState(object):
B=self.B,
clabel=self.clabel,
deg_corr=self.deg_corr,
max_BE=self.max_BE)
max_BE=self.max_BE,
ignore_degrees=self.ignore_degrees)
return state
def __setstate__(self, state):
......@@ -549,7 +555,9 @@ class BlockState(object):
S += libcommunity.deg_entropy_term(self.g._Graph__graph,
libcore.any(),
self.overlap_stats,
self.N)
self.N,
_prop("e", self.g, self.eweight),
_prop("v", self.g, self.ignore_degrees))
if multigraph:
S += libcommunity.entropy_parallel(self.g._Graph__graph,
......@@ -2181,6 +2189,7 @@ def minimize_blockmodel_dl(g, deg_corr=True, overlap=False, ec=None,
nested_overlap = kwargs.get("nested_overlap", False)
nonoverlap_compare = kwargs.get("nonoverlap_compare", False)
dl_ent = kwargs.get("dl_ent", False)
ignore_degrees = kwargs.get("ignore_degrees", None)
if minimize_state is None:
minimize_state = MinimizeState()
......@@ -2289,7 +2298,7 @@ def minimize_blockmodel_dl(g, deg_corr=True, overlap=False, ec=None,
else:
state = BlockState(g, B=g.num_vertices(), deg_corr=deg_corr,
vweight=vweight, eweight=eweight, clabel=clabel,
max_BE=max_BE)
max_BE=max_BE, ignore_degrees=ignore_degrees)
else:
if overlap:
......
src/graph_tool/community/nested_blockmodel.py
View file @ 36514eb1
......@@ -78,7 +78,7 @@ class NestedBlockState(object):
def __init__(self, g, eweight=None, vweight=None, ec=None, bs=None, Bs=None,
deg_corr=True, overlap=False, layers=False, clabel=None,
max_BE=1000):
max_BE=1000, **kwargs):
L = len(Bs) if Bs is not None else len(bs)
self.g = cg = g
self.vweight = vcount = vweight
......@@ -90,6 +90,7 @@ class NestedBlockState(object):
self.overlap = overlap
self.deg_corr = deg_corr
self.clabel = clabel if clabel is not None else g.new_vertex_property("int")
self.ignore_degrees = kwargs.get("ignore_degrees", None)
for l in range(L):
Bl = Bs[l] if Bs is not None else None
......@@ -117,7 +118,8 @@ class NestedBlockState(object):
vweight=vcount,
deg_corr=deg_corr != False,
#clabel=self.clabel,
max_BE=max_BE)
max_BE=max_BE,
ignore_degrees=self.ignore_degrees)
else:
state = CovariateBlockState(g, B=Bl, b=bl,
ec=ec,
......@@ -742,7 +744,8 @@ def replace_level(l, state, min_B=None, max_B=None, max_b=None, nsweeps=10,
##exaustive=g.num_vertices() <= 100,
#minimize_state=minimize_state.minimize_state, >>>>>> HERE <<<<<
checkpoint=checkpoint,
dl_ent=dl_ent)
dl_ent=dl_ent,
ignore_degrees=state.ignore_degrees if l == 0 else None)
if _bm_test():
assert (res.clabel.a == cclabel.a).all(), (res.clabel.a, cclabel.a)
......@@ -1240,6 +1243,7 @@ def init_nested_state(g, Bs, ec=None, deg_corr=True, overlap=False,
"""
dl_ent = kwargs.get("dl_ent", False)
ignore_degrees = kwargs.get("ignore_degrees", None)
if minimize_state is None:
minimize_state = NestedMinimizeState()
......@@ -1247,7 +1251,8 @@ def init_nested_state(g, Bs, ec=None, deg_corr=True, overlap=False,
state = NestedBlockState(g, ec=ec, layers=layers, eweight=eweight,
vweight=vweight, Bs=[1], deg_corr=deg_corr,
overlap=overlap, clabel=clabel)
overlap=overlap, clabel=clabel,
ignore_degrees=ignore_degrees)
chkp = get_checkpoint_wrap(checkpoint, state, minimize_state, dl_ent)
......@@ -1274,7 +1279,8 @@ def init_nested_state(g, Bs, ec=None, deg_corr=True, overlap=False,
eweight=ecount,
deg_corr=deg_corr != False,
#clabel=clabel,
max_BE=max_BE)
max_BE=max_BE,
ignore_degrees=ignore_degrees)
else:
if overlap:
if confine_layers:
......@@ -1570,6 +1576,7 @@ def minimize_nested_blockmodel_dl(g, Bs=None, bs=None, min_B=None,
"""
dl_ent = kwargs.get("dl_ent", False)
ignore_degrees = kwargs.get("ignore_degrees", None)
if minimize_state is None:
minimize_state = NestedMinimizeState()
......@@ -1659,12 +1666,13 @@ def minimize_nested_blockmodel_dl(g, Bs=None, bs=None, min_B=None,
sequential=sequential,
parallel=parallel,
checkpoint=checkpoint,
minimize_state=minimize_state, dl_ent=dl_ent)
minimize_state=minimize_state, dl_ent=dl_ent,
ignore_degrees=ignore_degrees)
else:
state = NestedBlockState(g, ec=ec, layers=layers, bs=bs,
deg_corr=deg_corr, overlap=overlap,
eweight=eweight, vweight=vweight,
clabel=clabel)
clabel=clabel, ignore_degrees=ignore_degrees)
minimize_state.sync(state)
......
src/graph_tool/community/overlap_blockmodel.py
View file @ 36514eb1
......@@ -581,8 +581,8 @@ class OverlapBlockState(BlockState):
if self.deg_corr:
S += libcommunity.deg_entropy_term(self.g._Graph__graph,
_prop("v", self.g, self.b),
self.overlap_stats,
self.N)
self.overlap_stats, self.N,
libcore.any(), libcore.any())
if self.deg_corr:
S -= E
else:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment