Commit 3172e99e authored by Tiago Peixoto's avatar Tiago Peixoto

uncertain_blockmodel.py: Fix issue with entropy computation

This also adds the "density" entropy argument.
parent f1193f1b
......@@ -180,20 +180,23 @@ struct Measured
return S;
}
double entropy()
double entropy(bool latent_edges, bool density)
{
double S = 0;
size_t gE = 0;
for (auto e : edges_range(_g))
if (latent_edges)
{
S += lbinom(_n[e], _x[e]);
gE++;
}
S += (_NP - gE) * lbinom(_n_default, _x_default);
for (auto e : edges_range(_g))
{
S += lbinom(_n[e], _x[e]);
gE++;
}
S += (_NP - gE) * lbinom(_n_default, _x_default);
S += get_MP(_T, _M);
S += get_MP(_T, _M);
}
if (_E_prior)
if (density && _E_prior)
S += _E * _pe - lgamma_fast(_E + 1) - exp(_pe);
return -S;
......@@ -213,13 +216,14 @@ struct Measured
double dS = _block_state.template modify_edge_dS<false>(source(e, _u),
target(e, _u),
e, _recs, ea);
if (ea.density && _E_prior)
{
dS += _pe;
dS += lgamma_fast(_E) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if (_E_prior)
{
dS += _pe;
dS += lgamma_fast(_E) - lgamma_fast(_E + 1);
}
if (_eweight[e] == 1 && (_self_loops || u != v))
{
auto& m = get_edge<false>(u, v);
......@@ -237,13 +241,14 @@ struct Measured
auto& e = get_u_edge(u, v);
double dS = _block_state.template modify_edge_dS<true>(u, v, e,
_recs, ea);
if (ea.density && _E_prior)
{
dS -= _pe;
dS += lgamma_fast(_E + 2) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if (_E_prior)
{
dS -= _pe;
dS += lgamma_fast(_E + 2) - lgamma_fast(_E + 1);
}
if ((e == _null_edge || _eweight[e] == 0) && (_self_loops || u != v))
{
auto& m = get_edge<false>(u, v);
......
......@@ -63,7 +63,8 @@ void export_uncertain_state()
class_<uentropy_args_t, bases<entropy_args_t>>("uentropy_args",
init<entropy_args_t>())
.def_readwrite("latent_edges", &uentropy_args_t::latent_edges);
.def_readwrite("latent_edges", &uentropy_args_t::latent_edges)
.def_readwrite("density", &uentropy_args_t::density);
def("make_uncertain_state", &make_uncertain_state);
......
......@@ -37,6 +37,7 @@ struct uentropy_args_t:
uentropy_args_t(const entropy_args_t& ea)
: entropy_args_t(ea){}
bool latent_edges;
bool density;
};
......@@ -143,34 +144,37 @@ struct Uncertain
return _get_edge<insert>(u, v, _g, _edges);
}
double entropy()
double entropy(bool latent_edges, bool density)
{
double S = 0;
for (auto m : edges_range(_g))
if (latent_edges)
{
double q_e = _q[m];
if (q_e == std::numeric_limits<double>::infinity())
continue;
auto& e = get_u_edge<false>(source(m, _g), target(m, _g));
if (e != _null_edge && _eweight[e] > 0 &&
(_self_loops || (source(e, _u) != target(e, _u))))
S += q_e;
}
for (auto m : edges_range(_g))
{
double q_e = _q[m];
if (q_e == std::numeric_limits<double>::infinity())
continue;
auto& e = get_u_edge<false>(source(m, _g), target(m, _g));
if (e != _null_edge && _eweight[e] > 0 &&
(_self_loops || (source(e, _u) != target(e, _u))))
S += q_e;
}
for (auto e : edges_range(_u))
{
auto& m = get_edge<false>(source(e, _u), target(e, _u));
if (m != _null_edge || _eweight[e] == 0 ||
(!_self_loops && source(m, _g) == target(m, _g)))
continue;
if (_q_default == std::numeric_limits<double>::infinity())
continue;
S += _q_default;
for (auto e : edges_range(_u))
{
auto& m = get_edge<false>(source(e, _u), target(e, _u));
if (m != _null_edge || _eweight[e] == 0 ||
(!_self_loops && source(m, _g) == target(m, _g)))
continue;
if (_q_default == std::numeric_limits<double>::infinity())
continue;
S += _q_default;
}
S += _S_const;
}
if (_E_prior)
if (density && _E_prior)
S += _E * _pe - lgamma_fast(_E + 1) - exp(_pe);
S += _S_const;
return -S;
}
......@@ -181,13 +185,14 @@ struct Uncertain
double dS = _block_state.template modify_edge_dS<false>(source(e, _u),
target(e, _u),
e, _recs, ea);
if (ea.density && _E_prior)
{
dS += _pe;
dS += lgamma_fast(_E) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if (_E_prior)
{
dS += _pe;
dS += lgamma_fast(_E) - lgamma_fast(_E + 1);
}
if (_eweight[e] == 1 && (_self_loops || u != v))
{
auto& m = get_edge<false>(u, v);
......@@ -203,13 +208,14 @@ struct Uncertain
auto& e = get_u_edge(u, v);
double dS = _block_state.template modify_edge_dS<true>(u, v, e,
_recs, ea);
if (ea.density && _E_prior)
{
dS -= _pe;
dS += lgamma_fast(_E + 2) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if (_E_prior)
{
dS -= _pe;
dS += lgamma_fast(_E + 2) - lgamma_fast(_E + 1);
}
if ((e == _null_edge || _eweight[e] == 0) && (_self_loops || u != v))
{
auto& m = get_edge<false>(u, v);
......
......@@ -34,9 +34,10 @@ from . nested_blockmodel import *
from . blockmodel import _bm_test
def get_uentropy_args(kargs):
ea = get_entropy_args(kargs, ignore=["latent_edges"])
ea = get_entropy_args(kargs, ignore=["latent_edges", "density"])
uea = libinference.uentropy_args(ea)
uea.latent_edges = kargs.get("latent_edges", True)
uea.density = kargs.get("density", True)
return uea
class UncertainBaseState(object):
......@@ -107,17 +108,37 @@ class UncertainBaseState(object):
self.slist = Vector_size_t(init=edges[:,0])
self.tlist = Vector_size_t(init=edges[:,1])
init_q_cache()
def get_block_state(self):
return self.bstate
def entropy(self, **kwargs):
def entropy(self, latent_edges=True, density=True, **kwargs):
if self.nbstate is None:
return self._state.entropy() + self.bstate.entropy(**kwargs)
S = self._state.entropy(latent_edges, density) + \
self.bstate.entropy(**kwargs)
else:
return self._state.entropy() + self.nbstate.entropy(**kwargs)
S = self._state.entropy(latent_edges, density) + \
self.nbstate.entropy(**kwargs)
if kwargs.get("test", True) and _bm_test():
assert not isnan(S) and not isinf(S), \
"invalid entropy %g (%s) " % (S, str(args))
args = kwargs.copy()
args["test"] = False
state_copy = self.copy()
Salt = state_copy.entropy(latent_edges, density, **args)
assert math.isclose(S, Salt, abs_tol=1e-8), \
"entropy discrepancy after copying (%g %g %g)" % (S, Salt,
S - Salt)
return S
def _algo_sweep(self, algo, r=.5, **kwargs):
kwargs = kwargs.copy()
beta = kwargs.get("beta", 1.)
niter = kwargs.get("niter", 1)
verbose = kwargs.get("verbose", False)
......@@ -126,32 +147,33 @@ class UncertainBaseState(object):
dentropy_args = dict(self.bstate._entropy_args,
**kwargs.get("entropy_args", {}))
entropy_args = get_uentropy_args(dentropy_args)
kwargs.get("entropy_args", {}).pop("latent_edges", None)
kwargs.get("entropy_args", {}).pop("density", None)
state = self._state
mcmc_state = DictState(locals())
if _bm_test():
Si = self.entropy(**dentropy_args)
try:
if self.nbstate is None:
self.bstate._clear_egroups()
else:
self.nbstate._clear_egroups()
if numpy.random.random() < r:
edges = True
dS, nattempts, nmoves = self._mcmc_sweep(mcmc_state)
else:
edges = False
if self.nbstate is None:
self.bstate._clear_egroups()
else:
self.nbstate._clear_egroups()
if numpy.random.random() < r:
edges = True
dS, nattempts, nmoves = self._mcmc_sweep(mcmc_state)
dS, nattempts, nmoves = algo(self.bstate, **kwargs)
else:
edges = False
if self.nbstate is None:
dS, nattempts, nmoves = algo(self.bstate, **kwargs)
else:
dS, nattempts, nmoves = algo(self.nbstate, **kwargs)
finally:
if _bm_test():
Sf = self.entropy(**dentropy_args)
assert math.isclose(dS, (Sf - Si), abs_tol=1e-8), \
"inconsistent entropy delta %g (%g): %s %s" % (dS, Sf - Si, edges,
str(dentropy_args))
dS, nattempts, nmoves = algo(self.nbstate, **kwargs)
if _bm_test():
Sf = self.entropy(**dentropy_args)
assert math.isclose(dS, (Sf - Si), abs_tol=1e-8), \
"inconsistent entropy delta %g (%g): %s %s" % (dS, Sf - Si, edges,
str(dentropy_args))
return dS, nattempts, nmoves
......
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