Skip to content
GitLab
Commit 500878c3 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Fix bug with waiting times

parent 104e09ac
Hide whitespace changes
Inline Side-by-side
src/graph/inference/graph_blockmodel.hh
View file @ 500878c3
......@@ -56,9 +56,9 @@ typedef mpl::vector2<simple_degs_t, degs_map_t> degs_tr;
((pclabel,, vmap_t, 0)) \
((merge_map,, vmap_t, 0)) \
((deg_corr,, bool, 0)) \
((use_waiting,, bool, 0)) \
((vtfield,, vprop_map_t<double>::type, 0)) \
((btfield,, vprop_map_t<double>::type, 0)) \
((ecov_type,, int, 0)) \
((ecov,, eprop_map_t<double>::type, 0)) \
((bcovsum,, vprop_map_t<double>::type, 0)) \
((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))
GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)
......@@ -97,8 +97,10 @@ public:
}
// remove a vertex from its current block
template <class EFilt>
void remove_vertex(size_t v, EFilt&& efilt)
template <class EFilt, class EOP>
void remove_vertex(size_t v,
EFilt&& efilt = [](auto&){ return false; },
EOP&& eop = [](auto&, auto&){})
{
typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
......@@ -132,6 +134,8 @@ public:
if (_mrs[me] == 0)
_emat.remove_me(r, s, me, _bg);
}
eop(e, me);
}
if (self_weight > 0)
......@@ -166,6 +170,8 @@ public:
if (_mrs[me] == 0)
_emat.remove_me(s, r, me, _bg);
eop(e, me);
}
_wr[r] -= _vweight[v];
......@@ -176,11 +182,20 @@ public:
if (is_partition_stats_enabled())
get_partition_stats(v).remove_vertex(v, r, _deg_corr, _g, _vweight,
_eweight, _degs);
if (_ecov_type == 3) // waiting times
{
if (_ignore_degrees[v] > 0)
{
double dt = out_degreeS()(v, _g, _ecov);
_bcovsum[r] -= dt;
}
}
}
void remove_vertex(size_t v)
{
remove_vertex(v, [](auto&){ return false; });
remove_vertex(v, [](auto&){ return false; }, [](auto&, auto&){});
}
template <class Vlist>
......@@ -203,7 +218,8 @@ public:
}
for (auto v : vset)
remove_vertex(v, [&](auto& e) { return eset.find(e) != eset.end(); });
remove_vertex(v, [&](auto& e) { return eset.find(e) != eset.end(); },
[](auto&, auto&){});
for (auto& e : eset)
{
......@@ -234,8 +250,10 @@ public:
}
// add a vertex to block r
template <class BEdge, class Efilt>
void add_vertex(size_t v, size_t r, BEdge&& bedge, Efilt&& efilt)
template <class BEdge, class Efilt, class EOP>
void add_vertex(size_t v, size_t r, BEdge&& bedge,
Efilt&& efilt = [](auto&){ return false; },
EOP&& eop = [](auto&, auto&){})
{
typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
typedef typename graph_traits<bg_t>::edge_descriptor bedge_t;
......@@ -279,6 +297,8 @@ public:
_mrp[r] += ew;
_mrm[s] += ew;
}
eop(e, me);
}
if (self_weight > 0)
......@@ -322,6 +342,8 @@ public:
_mrp[s] += ew;
_mrm[r] += ew;
eop(e, me);
}
_wr[r] += _vweight[v];
......@@ -333,11 +355,23 @@ public:
if (is_partition_stats_enabled())
get_partition_stats(v).add_vertex(v, r, _deg_corr, _g, _vweight,
_eweight, _degs);
if (_ecov_type == 3) // waiting times
{
if (_ignore_degrees[v] > 0)
{
double dt = out_degreeS()(v, _g, _ecov);
_bcovsum[r] += dt;
}
}
}
void add_vertex(size_t v, size_t r)
{
add_vertex(v, r, _emat.get_bedge_map(), [](auto&) { return false; });
add_vertex(v, r, _emat.get_bedge_map(),
[](auto&){ return false; },
[](auto&, auto&){});
}
template <class Vlist, class Blist>
......@@ -368,7 +402,8 @@ public:
for (auto vr : vset)
add_vertex(vr.first, vr.second, bedge,
[&](auto& e){ return eset.find(e) != eset.end(); });
[&](auto& e){ return eset.find(e) != eset.end(); },
[](auto&, auto&){});
for (auto e : eset)
{
......@@ -798,20 +833,23 @@ public:
dS += ps.get_delta_edges_dl(v, r, nr, _vweight, _g);
}
if (_use_waiting && (r != nr) && _ignore_degrees[v] > 0)
if (_ecov_type == 3) // waiting times
{
double dt = _vtfield[v];
int k = out_degreeS()(v, _g, _eweight);
if (_mrp[r] > 0)
dS -= _mrp[r] * log(_btfield[r]) - lgamma_fast(_mrp[r]);
if (_mrp[nr] > 0)
dS -= _mrp[nr] * log(_btfield[nr]) - lgamma_fast(_mrp[nr]);
if (_mrp[r] > k)
dS += (_mrp[r] - k) * log(_btfield[r] - dt)
- lgamma_fast(_mrp[r] - k);
if (_mrp[nr] + k > 0)
dS += (_mrp[nr] + k) * log(_btfield[nr] + dt)
- lgamma_fast(_mrp[nr] + k);
if ((r != nr) && _ignore_degrees[v] > 0)
{
double dt = out_degreeS()(v, _g, _ecov);
int k = out_degreeS()(v, _g, _eweight);
if (_mrp[r] > 0)
dS -= _mrp[r] * log(_bcovsum[r]) - lgamma_fast(_mrp[r]);
if (_mrp[nr] > 0)
dS -= _mrp[nr] * log(_bcovsum[nr]) - lgamma_fast(_mrp[nr]);
if (_mrp[r] > k)
dS += (_mrp[r] - k) * log(_bcovsum[r] - dt)
- lgamma_fast(_mrp[r] - k);
if (_mrp[nr] + k > 0)
dS += (_mrp[nr] + k) * log(_bcovsum[nr] + dt)
- lgamma_fast(_mrp[nr] + k);
}
}
return dS;
......@@ -1078,11 +1116,11 @@ public:
else
S = sparse_entropy(multigraph, deg_entropy);
if (_use_waiting)
if (_ecov_type == 3) // waiting times
{
for (auto r : vertices_range(_bg))
if (_btfield[r] > 0)
S += _mrp[r] * log(_btfield[r]) - lgamma_fast(_mrp[r]);
if (_bcovsum[r] > 0)
S += _mrp[r] * log(_bcovsum[r]) - lgamma_fast(_mrp[r]);
}
return S;
}
......
src/graph_tool/inference/blockmodel.py
View file @ 500878c3
......@@ -49,14 +49,14 @@ def _bm_test():
global __test__
return __test__
def get_block_graph(g, B, b, vcount=None, ecount=None, edt=None):
def get_block_graph(g, B, b, vcount=None, ecount=None, ecov=None):
if isinstance(ecount, libinference.unity_eprop_t):
ecount = None
if isinstance(vcount, libinference.unity_vprop_t):
vcount = None
aeprops = []
if edt is not None:
aeprops.append(edt)
if ecov is not None:
aeprops.append(ecov)
cg, br, vcount, ecount, av, ae = condensation_graph(g, b,
vweight=vcount,
eweight=ecount,
......@@ -64,8 +64,8 @@ def get_block_graph(g, B, b, vcount=None, ecount=None, edt=None):
self_loops=True)
cg.vp.count = vcount
cg.ep.count = ecount
if edt is not None:
cg.ep.dt = ae[0]
if ecov is not None:
cg.ep.ecov = ae[0]
del ae[0]
cg = Graph(cg, vorder=br)
......@@ -105,9 +105,9 @@ class BlockState(object):
the block graph. Otherwise a dense matrix will be used.
"""
def __init__(self, g, eweight=None, vweight=None, b=None, B=None,
clabel=None, pclabel=None, deg_corr=True, max_BE=1000,
**kwargs):
def __init__(self, g, eweight=None, vweight=None, b=None, B=None, ecov=None,
ecov_type=None, clabel=None, pclabel=None, deg_corr=True,
max_BE=1000, **kwargs):
kwargs = kwargs.copy()
# initialize weights to unity, if necessary
......@@ -170,7 +170,7 @@ class BlockState(object):
# Construct block-graph
self.bg = get_block_graph(g, B, self.b, self.vweight, self.eweight,
edt=kwargs.get("edt", None))
ecov=ecov)
self.bg.set_fast_edge_removal()
self.mrs = self.bg.ep["count"]
......@@ -232,18 +232,17 @@ class BlockState(object):
self.block_list = Vector_size_t()
self.block_list.extend(arange(self.B, dtype="int"))
self.edt = extract_arg(kwargs, "edt", None)
self.use_waiting = self.edt is not None
if self.use_waiting:
self.vtfield = self.g.degree_property_map("out", self.edt)
self.btfield = self.bg.degree_property_map("out", self.bg.ep.dt)
self.ecov = ecov
self.ecov_type = 3 if self.ecov is not None else 0
if self.ecov_type == 3: # waiting times
self.bcovsum = self.bg.degree_property_map("out", self.bg.ep.ecov)
tokens = self.ignore_degrees.fa == 0
token_groups = bincount(self.b.fa[tokens]) > 0
token_groups.resize(self.bg.num_vertices())
self.btfield.a[token_groups] = 0
self.bcovsum.a[token_groups] = 0
else:
self.vtfield = self.g.new_vp("double")
self.btfield = self.bg.new_vp("double")
self.ecov = self.g.new_ep("double")
self.bcovsum = self.bg.new_vp("double")
self._abg = self.bg._get_any()
self._state = libinference.make_block_state(self, _get_rng())
......@@ -284,7 +283,7 @@ class BlockState(object):
ignore_degrees=kwargs.pop("ignore_degrees", self.ignore_degrees),
degs=self.degs.copy(),
merge_map=kwargs.get("merge_map", self.merge_map.copy()),
edt=self.edt,
ecov=self.ecov,
**dmask(kwargs, ["ignore_degrees", "merge_map"]))
else:
state = OverlapBlockState(self.g if g is None else g,
......
Supports Markdown
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