Commit db70b049 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Fix minor issues with edge prediction and model averaging

parent 3cd04a82
......@@ -788,6 +788,7 @@ public:
for (auto e : edges_range(bg))
{
assert(get_me(source(e, bg),target(e, bg)) == _null_edge);
_mat[source(e, bg)][target(e, bg)] = e;
if (!is_directed::apply<BGraph>::type::value)
_mat[target(e, bg)][source(e, bg)] = e;
......@@ -895,7 +896,10 @@ public:
_hash.resize(num_vertices(bg), ehash_t(0, _hash_function));
for (auto e : edges_range(bg))
{
assert(get_me(source(e, bg), target(e, bg)) == _null_edge);
put_me(source(e, bg), target(e, bg), e);
}
auto bedge_c = _bedge.get_checked();
for (auto e : edges_range(g))
......
......@@ -277,9 +277,6 @@ class BlockState(object):
self.mrs = self.bg.ep["count"]
self.wr = self.bg.vp["count"]
del self.bg.ep["count"]
del self.bg.vp["count"]
self.mrp = self.bg.degree_property_map("out", weight=self.mrs)
if g.is_directed():
......@@ -970,11 +967,12 @@ class BlockState(object):
return self._state.get_move_prob(int(v), s, self.b[v], c, True)
def get_edges_prob(self, edge_list, missing=True, entropy_args={}):
"""Compute the log-probability of the missing (or spurious if ``missing=False``)
edges given by ``edge_list`` (a list of ``(source, target)`` tuples, or
:meth:`~graph_tool.Edge` instances). The values in ``entropy_args`` are
passed to :meth:`graph_tool.BlockState.entropy()` to calculate the
log-probability.
"""Compute the unnormalized log-probability of the missing (or spurious if
``missing=False``) edges given by ``edge_list`` (a list of ``(source,
target)`` tuples, or :meth:`~graph_tool.Edge` instances). The values in
``entropy_args`` are passed to :meth:`graph_tool.BlockState.entropy()`
to calculate the log-probability.
"""
pos = {}
for u, v in edge_list:
......@@ -989,9 +987,14 @@ class BlockState(object):
if missing:
new_es = []
for u, v in edge_list:
e = self.g.add_edge(u, v)
if self.is_weighted:
self.eweight[e] = 1
if not self.is_weighted:
e = self.g.add_edge(u, v)
else:
e = self.g.edge(u, v)
if e is None:
e = self.g.add_edge(u, v)
self.eweight[e] = 0
self.eweight[e] += 1
new_es.append(e)
else:
old_es = []
......@@ -1019,8 +1022,14 @@ class BlockState(object):
finally:
if missing:
for e in new_es:
self.g.remove_edge(e)
if self.is_weighted:
for e in reversed(new_es):
self.eweight[e] -= 1
if self.eweight[e] == 0:
self.g.remove_edge(e)
else:
for e in reversed(new_es):
self.g.remove_edge(e)
else:
for u, v in old_es:
if self.is_weighted:
......@@ -1031,12 +1040,22 @@ class BlockState(object):
self.eweight[e] += 1
else:
self.g.add_edge(u, v)
self.add_vertex(pos.keys(), pos.values())
if missing:
return Si - Sf
else:
return Sf - Si
L = Si - Sf
if _bm_test():
state = self.copy()
set_test(False)
L_alt = state.get_edges_prob(edge_list, missing=missing,
entropy_args=entropy_args)
set_test(True)
assert abs(L - L_alt) < 1e-6, \
"inconsistent missing=%s edge probability (%g, %g): %s, %s" % \
(str(missing), L, L_alt, str(entropy_args), str(edge_list))
return L
def _mcmc_sweep_dispatch(self, mcmc_state):
return libinference.mcmc_sweep(mcmc_state, self._state,
......@@ -1582,7 +1601,7 @@ class BlockState(object):
Returns
-------
p : :class:`~graph_tool.libcore.any`
p : :class:`~graph_tool.PropertyMap`
Edge property map with updated edge marginals.
Examples
......@@ -1860,13 +1879,13 @@ def bethe_entropy(g, p):
.. math::
H = -\sum_{e,(r,s)}\pi_{(r,s)}^e\ln\pi_{(r,s)}^e - \sum_{v,r}(1-k_i)\pi_r^v\ln\pi_r^v,
H = -\sum_{ij}A_{ij}\sum_{rs}\pi_{ij}(r,s)\ln\pi_{ij}(r,s) - \sum_i(1-k_i)\sum_r\pi_i(r)\ln\pi_i(r),
where :math:`\pi_{(r,s)}^e` is the marginal probability that the endpoints
of the edge :math:`e` belong to blocks :math:`(r,s)`, and :math:`\pi_r^v` is
the marginal probability that vertex :math:`v` belongs to block :math:`r`,
and :math:`k_i` is the degree of vertex :math:`v` (or total degree for
directed graphs).
where :math:`\pi_{ij}(r,s)` is the marginal probability that vertices
:math:`i` and :math:`j` belong to blocks :math:`r` and :math:`s`,
respectively, and :math:`\pi_i(r)` is the marginal probability that vertex
:math:`i` belongs to block :math:`r`, and :math:`k_i` is the degree of
vertex :math:`i` (or total degree for directed graphs).
References
----------
......@@ -1906,9 +1925,9 @@ def mf_entropy(g, p):
.. math::
H = - \sum_{v,r}\pi_r^v\ln\pi_r^v,
H = - \sum_{i,r}\pi_i(r)ln\pi_i(r),
where :math:`\pi_r^v` is the marginal probability that vertex :math:`v`
where :math:`\pi_i(r)` is the marginal probability that vertex :math:`i`
belongs to block :math:`r`.
References
......
......@@ -477,7 +477,7 @@ class LayeredBlockState(OverlapBlockState, BlockState):
return bg, mrs, ec, rec, drec
def get_block_state(self, b=None, vweight=False, deg_corr=False,
overlap=False, layers=None, **kwargs):
overlap=False, layers=True, **kwargs):
r"""Returns a :class:`~graph_tool.inference.LayeredBlockState`` corresponding
to the block graph. The parameters have the same meaning as the in the
constructor."""
......@@ -749,10 +749,19 @@ class LayeredBlockState(OverlapBlockState, BlockState):
self.ec[e] = l[0]
self.add_vertex(pos.keys(), pos.values())
if missing:
return Si - Sf
else:
return Sf - Si
L = Si - Sf
if _bm_test():
state = self.copy()
set_test(False)
L_alt = state.get_edges_prob(edge_list, missing=missing,
entropy_args=entropy_args)
set_test(True)
assert abs(L - L_alt) < 1e-6, \
"inconsistent missing=%s edge probability (%g, %g): %s, %s" % \
(str(missing), L, L_alt, str(entropy_args), str(edge_list))
return L
def _mcmc_sweep_dispatch(self, mcmc_state):
if not self.overlap:
......
......@@ -28,7 +28,7 @@ from .. import Vector_size_t, Vector_double
import numpy
from . util import *
def mcmc_equilibrate(state, wait=10, nbreaks=2, max_niter=numpy.inf,
def mcmc_equilibrate(state, wait=1000, nbreaks=2, max_niter=numpy.inf,
force_niter=None, epsilon=0, gibbs=False,
block_moves=False, mcmc_args={}, entropy_args={},
history=False, callback=None, verbose=False):
......@@ -38,7 +38,7 @@ def mcmc_equilibrate(state, wait=10, nbreaks=2, max_niter=numpy.inf,
----------
state : Any state class (e.g. :class:`~graph_tool.inference.BlockState`)
Initial state. This state will be modified during the algorithm.
wait : ``int`` (optional, default: ``10``)
wait : ``int`` (optional, default: ``1000``)
Number of iterations to wait for a record-breaking event.
nbreaks : ``int`` (optional, default: ``2``)
Number of iteration intervals (of size ``wait``) without record-breaking
......
......@@ -66,6 +66,7 @@ class NestedBlockState(object):
**kwargs : keyword arguments
Keyword arguments to be passed to base type constructor.
"""
def __init__(self, g, bs, base_type=BlockState, hstate_args={},
hentropy_args={}, sampling=False, **kwargs):
self.g = g
......@@ -121,15 +122,16 @@ class NestedBlockState(object):
g = copy.deepcopy(self.g, memo)
return self.copy(g=g)
def get_bs(self):
return [s.b.fa for s in self.levels]
def copy(self, g=None, bs=None, **kwargs):
def copy(self, g=None, bs=None, hstate_args=None, hentropy_args=None,
sampling=None, **kwargs):
r"""Copies the block state. The parameters override the state properties,
and have the same meaning as in the constructor."""
bs = self.get_bs() if bs is None else bs
return NestedBlockState(self.g if g is None else g, bs,
base_type=type(self.levels[0]),
hstate_args=self.hstate_args if hstate_args is None else hstate_args,
hentropy_args=self.hentropy_args if hentropy_args is None else hentropy_args,
sampling=self.sampling if sampling is None else sampling,
**overlay(self.kwargs, **kwargs))
def __getstate__(self):
......@@ -143,6 +145,17 @@ class NestedBlockState(object):
conv_pickle_state(state)
self.__init__(**overlay(dmask(state, ["kwargs"]), **state["kwargs"]))
def get_bs(self):
"""Get hierarchy levels as a list of :class:`numpy.ndarray` objects with the
group memberships at each level.
"""
return [s.b.fa for s in self.levels]
def get_levels(self):
"""Get hierarchy levels as a list of :class:`~graph_tool.inference.BlockState`
instances."""
return self.levels
def project_partition(self, j, l):
"""Project partition of level ``j`` onto level ``l``, and return it."""
b = self.levels[l].b.copy()
......@@ -233,8 +246,7 @@ class NestedBlockState(object):
bstate = self.levels[l]
if l > 0:
eargs = overlay(self.hentropy_args,
**overlay(kwargs, multigraph=True))
eargs = self.hentropy_args
else:
eargs = kwargs
......@@ -295,16 +307,28 @@ class NestedBlockState(object):
log-probability.
"""
L = 0
for l, state in enumerate(self.levels):
eargs = overlay(entropy_args,
edges_dl=(l == (len(self.levels) - 1)))
for l, lstate in enumerate(self.levels):
if l > 0:
eargs = overlay(eargs, **self.hentropy_args)
L += state.get_edges_prob(edge_list, missing=missing,
entropy_args=eargs)
edge_list = [(state.b[e[0]], state.b[e[1]]) for e in (tuple(e_) for e_ in edge_list)]
eargs = self.hentropy_args
else:
eargs = entropy_args
eargs = overlay(eargs, dl=True,
edges_dl=(l == (len(self.levels) - 1)))
if self.sampling:
lstate._couple_state(None, None)
if l > 0:
lstate._state.sync_emat()
L += lstate.get_edges_prob(edge_list, missing, entropy_args=eargs)
if isinstance(self.levels[0], LayeredBlockState):
edge_list = [(lstate.b[u], lstate.b[v], l) for u, v, l in edge_list]
else:
edge_list = [(lstate.b[u], lstate.b[v]) for u, v in (tuple(e_) for e_ in edge_list)]
return L
def get_bstack(self):
"""Return the nested levels as individual graphs.
......@@ -491,6 +515,12 @@ class NestedBlockState(object):
The arguments accepted are the same as in
:method:`graph_tool.inference.BlockState.mcmc_sweep`.
"""
c = kwargs.get("c", None)
if c is None:
c = [1] + [numpy.inf] * (len(self.levels) - 1)
kwargs = kwargs.copy()
kwargs["c"] = c
return self._h_sweep(lambda s, **a: s.mcmc_sweep(**a), **kwargs)
def gibbs_sweep(self, **kwargs):
......@@ -511,27 +541,6 @@ class NestedBlockState(object):
"""
return self._h_sweep(lambda s, **a: s.multicanonical_sweep(**a))
def get_edges_prob(self, edge_list, missing=True, entropy_args={}):
"""Compute the log-probability of the missing (or spurious if ``missing=False``)
edges given by ``edge_list`` (a list of ``(source, target)`` tuples, or
:meth:`~graph_tool.Edge` instances). The values in ``entropy_args`` are
passed to :meth:`graph_tool.NestedBlockState.entropy()` to calculate the
log-probability.
"""
S = 0
for l, lstate in enumerate(self.levels):
if l > 0:
eargs = overlay(self.hentropy_args,
edges_dl=(l == len(self.levels) - 1))
else:
eargs = entropy_args
S += lstate.get_edges_prob(edge_list, missing, entropy_args=eargs)
if isinstance(self.levels[0], LayeredBlockState):
edge_list = [(lstate.b[u], lstate.b[v], l) for u, v, l in edge_list]
else:
edge_list = [(lstate.b[u], lstate.b[v]) for u, v in edge_list]
return S
def collect_partition_histogram(self, h=None, update=1):
r"""Collect a histogram of partitions.
......@@ -843,7 +852,6 @@ def get_hierarchy_tree(state, empty_branches=True):
t.set_vertex_filter(t.own_property(filt[0].copy()),
filt[1])
label = t.vertex_index.copy("int")
t.vp.mask = t.new_vp("bool", True)
order = t.own_property(g.degree_property_map("total").copy())
t_vertices = list(t.vertices())
......@@ -857,11 +865,6 @@ def get_hierarchy_tree(state, empty_branches=True):
t_vertices.append(t.add_vertex(s.num_vertices()))
label.a[-s.num_vertices():] = arange(s.num_vertices())
if "count" in s.vp:
t.vp.mask.a[-s.num_vertices():] = s.vp.count.fa > 0
else:
t.vp.mask.a[-s.num_vertices():] = True
# relative ordering based on total degree
count = s.ep["count"].copy("double")
for e in s.edges():
......@@ -879,7 +882,10 @@ def get_hierarchy_tree(state, empty_branches=True):
for vi, v in enumerate(g.vertices()):
w = t_vertices[vi + last_pos]
u = t_vertices[b[v] + pos]
if s.num_vertices() == 1:
u = t_vertices[pos]
else:
u = t_vertices[b[v] + pos]
t.add_edge(u, w)
last_pos = pos
......@@ -890,17 +896,17 @@ def get_hierarchy_tree(state, empty_branches=True):
if not empty_branches:
vmask = t.new_vertex_property("bool", True)
for vi in range(state.g.num_vertices(), t.num_vertices()):
v = t_vertices[vi]
t = GraphView(t, vfilt=vmask)
vmask = t.get_vertex_filter()[0]
for vi in range(state.g.num_vertices(ignore_filter=True),
t.num_vertices()):
v = t.vertex(t_vertices[vi])
if v.out_degree() == 0:
vmask[v] = False
while v.in_degree() > 0:
v = list(v.in_neighbours())[0]
vmask[v] = False
vmask[v] = True
t = GraphView(t, vfilt=logical_and(vmask.fa, t.vp.mask.fa))
t.vp.label = t.own_property(label)
t.vp.order = t.own_property(order)
t.vp.label = label
t.vp.order = order
t = Graph(t, prune=True)
label = t.vp.label
order = t.vp.order
......
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