Commit 31ae4357 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

inference.blockmodel: Slight change to signed edge covariate model

parent 03b4d2c0
This diff is collapsed.
......@@ -63,6 +63,7 @@ typedef mpl::vector2<std::true_type, std::false_type> use_hash_tr;
((wparams,, std::vector<std::vector<double>>, 0)) \
((recdx, &, std::vector<double>&, 0)) \
((Lrecdx, &, std::vector<double>&, 0)) \
((epsilon, &, std::vector<double>&, 0)) \
((allow_empty,, bool, 0))
GEN_STATE_BASE(OverlapBlockStateVirtualBase, OVERLAP_BLOCK_STATE_params)
......@@ -316,6 +317,9 @@ public:
}
else
{
for (size_t i = 0; i < this->_rec_types.size(); ++i)
_dBdx[i] = 0;
auto mid_op =
[&](auto& e, auto& me)
{
......@@ -334,7 +338,7 @@ public:
one * this->_rec[i][e]), 2) /
(mrs + one * this->_rec[0][e])));
_recdx[i] += dx;
_Lrecdx[i+1] += dx;
_dBdx[i] += dx;
}
}
......@@ -352,7 +356,6 @@ public:
(this->_bdrec[i][me] -
std::pow(this->_brec[i][me], 2) / mrs);
_recdx[i] -= dx;
_Lrecdx[i+1] -= dx;
}
}
......@@ -384,7 +387,13 @@ public:
}
};
for (size_t i = 0; i < _rec_types.size(); ++i)
_Lrecdx[i+1] -= _recdx[i] * _B_E_D;
modify_vertex<Add>(v, r, mid_op, coupled_end_op);
for (size_t i = 0; i < _rec_types.size(); ++i)
_Lrecdx[i+1] += _recdx[i] * _B_E_D;
}
}
}
......@@ -658,12 +667,14 @@ public:
positive_entries_op(i,
[&](auto N, auto x)
{ return positive_w_log_P(N, x, wp[0],
wp[1]);
wp[1],
this->_epsilon[i]);
},
[&](size_t B_E)
{ return positive_w_log_P(B_E,
_recsum[i],
wp[0], wp[1]);
wp[0], wp[1],
this->_epsilon[i]);
});
break;
case weight_type::DISCRETE_GEOMETRIC:
......@@ -727,12 +738,14 @@ public:
auto dx2 = get<2>(delta)[i];
dS -= -signed_w_log_P(ers, xrs, x2rs,
wp[0], wp[1],
wp[2], wp[3]);
wp[2], wp[3],
this->_epsilon[i]);
dS += -signed_w_log_P(ers + d,
xrs + dx,
x2rs + dx2,
wp[0], wp[1],
wp[2], wp[3]);
wp[2], wp[3],
this->_epsilon[i]);
if (std::isnan(wp[0]) &&
std::isnan(wp[1]))
{
......@@ -769,19 +782,22 @@ public:
{
dS -= -signed_w_log_P(_B_E, _recsum[i],
_recx2[i], wp[0], wp[1],
wp[2], wp[3]);
wp[2], wp[3], _epsilon[i]);
dS += -signed_w_log_P(_B_E + dB_E, _recsum[i],
_recx2[i] + dBx2, wp[0],
wp[1], wp[2], wp[3]);
wp[1], wp[2], wp[3],
_epsilon[i]);
}
if (dB_E_D != 0 || _dBdx[i] != 0)
{
dS -= -positive_w_log_P_alt(_B_E_D, _recdx[i],
wp[2], wp[3]);
dS += -positive_w_log_P_alt(_B_E_D + dB_E_D,
_recdx[i] + _dBdx[i],
wp[2], wp[3]);
dS -= -positive_w_log_P(_B_E_D, _recdx[i],
wp[2], wp[3],
_epsilon[i]);
dS += -positive_w_log_P(_B_E_D + dB_E_D,
_recdx[i] + _dBdx[i],
wp[2], wp[3],
_epsilon[i]);
}
if (dL == 0)
......@@ -792,14 +808,26 @@ public:
dL--;
}
if (_coupled_state == nullptr && _Lrecdx[0] >= 0)
if (_Lrecdx[0] >= 0)
{
size_t L = _Lrecdx[0];
dS -= -positive_w_log_P_alt(L, _Lrecdx[i+1],
wp[2], wp[3]);
dS += -positive_w_log_P_alt(L + dL,
_Lrecdx[i+1] + _dBdx[i],
wp[2], wp[3]);
size_t N_B_E_D = _B_E_D + dB_E_D;
dS -= -safelog(_B_E_D);
dS += -safelog(N_B_E_D);
_dBdx[i] = _recdx[i] * dB_E_D + _dBdx[i] * N_B_E_D;
if (_coupled_state == nullptr)
{
size_t L = _Lrecdx[0];
dS -= -positive_w_log_P(L, _Lrecdx[i+1],
wp[2], wp[3],
_epsilon[i]);
dS += -positive_w_log_P(L + dL,
_Lrecdx[i+1] +
_dBdx[i], wp[2],
wp[3], _epsilon[i]);
}
}
}
}
......@@ -1127,10 +1155,12 @@ public:
{
auto ers = _brec[0][me];
auto xrs = _brec[i][me];
S += -positive_w_log_P(ers, xrs, wp[0], wp[1]);
S += -positive_w_log_P(ers, xrs, wp[0], wp[1],
_epsilon[i]);
}
if (edges_dl && std::isnan(wp[0]) && std::isnan(wp[1]))
S += -positive_w_log_P(_B_E, _recsum[i], wp[0], wp[1]);
S += -positive_w_log_P(_B_E, _recsum[i], wp[0], wp[1],
_epsilon[i]);
break;
case weight_type::DISCRETE_GEOMETRIC:
for (auto me : edges_range(_bg))
......@@ -1173,15 +1203,16 @@ public:
auto xrs = _brec[i][me];
auto x2rs = _bdrec[i][me];
S += -signed_w_log_P(ers, xrs, x2rs, wp[0], wp[1],
wp[2], wp[3]);
wp[2], wp[3], _epsilon[i]);
}
if (std::isnan(wp[0]) && std::isnan(wp[1]))
{
if (edges_dl)
S += -signed_w_log_P(_B_E, _recsum[i], _recx2[i],
wp[0], wp[1], wp[2], wp[3]);
S += -positive_w_log_P_alt(_B_E_D, _recdx[i], wp[2],
wp[3]);
wp[0], wp[1], wp[2], wp[3],
_epsilon[i]);
S += -positive_w_log_P(_B_E_D, _recdx[i], wp[2], wp[3],
_epsilon[i]);
}
break;
case weight_type::DELTA_T: // waiting times
......
......@@ -280,13 +280,14 @@ inline double eterm_dense(size_t r, size_t s, int ers, double wr_r,
// exponential
template <class DT>
double positive_w_log_P(DT N, double x, double alpha, double beta)
double positive_w_log_P(DT N, double x, double alpha, double beta,
double epsilon)
{
if (N == 0)
return 0.;
if (std::isnan(alpha) && std::isnan(beta))
{
if (x == 0 || N == 1)
if (x < epsilon || N == 1)
return 0.;
else
return lgamma(N) - N * log(x);
......@@ -295,25 +296,17 @@ double positive_w_log_P(DT N, double x, double alpha, double beta)
(alpha + N) * log(beta + x);
}
template <class DT>
double positive_w_log_P_alt(DT N, double x, double alpha, double beta)
{
if (abs(x) < 1e-10)
x = 0;
return positive_w_log_P(N, x, alpha, beta);
}
// normal
template <class DT>
double signed_w_log_P(DT N, double x, double x2, double m0, double k0, double v0,
double nu0)
double nu0, double epsilon)
{
if (N == 0)
return 0.;
if (std::isnan(m0) && std::isnan(k0))
{
auto smu1 = x * (x / N);
if (N <= 2 || smu1 >= x2 || boost::math::relative_difference(x2, smu1) < 1e-10)
if (N <= 2 || smu1 >= x2 || (x2 - smu1) < std::pow(epsilon, 2))
return 0.;
else
return (lgamma((N - 1) / 2.) + log(x2) / 2. - ((N - 2) / 2.) *
......
......@@ -431,6 +431,11 @@ class BlockState(object):
self.Lrecdx = libcore.Vector_double(len(self.rec)+1)
self.Lrecdx[0] = -1
self.Lrecdx.resize(len(self.rec)+1)
self.epsilon = kwargs.pop("epsilon", libcore.Vector_double(len(self.rec)))
for i in range(len(self.rec)):
idx = self.rec[i].a != 0
if numpy.any(idx):
self.epsilon[i] = abs(self.rec[i].a[idx]).min() / 10
self.allow_empty = allow_empty
self._abg = self.bg._get_any()
......@@ -611,6 +616,7 @@ class BlockState(object):
allow_empty=kwargs.pop("allow_empty",
self.allow_empty),
Lrecdx=kwargs.pop("Lrecdx", self.Lrecdx.copy()),
epsilon=self.epsilon.copy(),
**kwargs)
else:
state = OverlapBlockState(self.g if g is None else g,
......@@ -629,6 +635,7 @@ class BlockState(object):
self.allow_empty),
max_BE=self.max_BE if max_BE is None else max_BE,
Lrecdx=kwargs.pop("Lrecdx", self.Lrecdx.copy()),
epsilon=self.epsilon.copy(),
**kwargs)
if self._coupled_state is not None:
......@@ -773,6 +780,7 @@ class BlockState(object):
clabel=kwargs.pop("clabel", self.get_bclabel()),
pclabel=kwargs.pop("pclabel", self.get_bpclabel()),
max_BE=self.max_BE,
epsilon=self.epsilon.copy(),
**kwargs)
if copy_coupled and self._coupled_state is not None:
......
......@@ -187,6 +187,8 @@ class LayeredBlockState(OverlapBlockState, BlockState):
self.drec = agg_state.drec
self.rec_types = agg_state.rec_types
self.rec_params = agg_state.rec_params
self.Lrecdx = agg_state.Lrecdx
self.epsilon = agg_state.epsilon
self.b = agg_state.b
self.B = agg_state.B
......
......@@ -33,12 +33,6 @@ from numpy import *
import numpy
import copy
def get_edges_dl(state, hstate_args, hentropy_args):
bclabel = state.get_bclabel()
bstate = state.get_block_state(b=bclabel, **dict(hstate_args,
Lrecdx=state.Lrecdx))
return bstate.entropy(**hentropy_args)
class NestedBlockState(object):
r"""The nested stochastic block model state of a given graph.
......@@ -114,7 +108,8 @@ class NestedBlockState(object):
self.Lrecdx.a = 0
self.Lrecdx[0] = len([s for s in self.levels if s._state.get_B_E_D() > 0])
for s in self.levels:
self.Lrecdx.a[1:] += s.recdx.a
self.Lrecdx.a[1:] += s.recdx.a * s._state.get_B_E_D()
s.epsilon.a = self.levels[0].epsilon.a
def _regen_levels(self):
for l in range(1, len(self.levels)):
......@@ -278,12 +273,24 @@ class NestedBlockState(object):
return S
def _Lrecdx_entropy(self, Lrecdx=None):
S_D = 0
if Lrecdx is None:
Lrecdx = self.Lrecdx
for s in self.levels:
B_E_D = s._state.get_B_E_D()
if B_E_D > 0:
S_D -= log(B_E_D)
S = 0
for i in range(len(Lrecdx)-1):
for i in range(len(self.levels[0].rec)):
if self.levels[0].rec_types[i] != libinference.rec_type.real_normal:
continue
S += -libinference.positive_w_log_P(Lrecdx[0], Lrecdx[i+1],
numpy.nan, numpy.nan)
numpy.nan, numpy.nan,
self.levels[0].epsilon[i])
S += S_D
return S
def entropy(self, **kwargs):
......@@ -305,7 +312,7 @@ class NestedBlockState(object):
Salt = state.entropy(test=False, **kwargs)
assert math.isclose(S, Salt, abs_tol=1e-8), \
"inconsistent entropy after copying (%g, %g, %g): %s" % \
(S, Salt, abs(S-Salt), str(kwargs))
(S, Salt, S-Salt, str(kwargs))
return S
......@@ -454,11 +461,27 @@ class NestedBlockState(object):
def callback(s):
S = 0
if l < len(self.levels) - 1:
S += get_edges_dl(s, self.hstate_args,
dict(self.hentropy_args,
edges_dl=(l + 1 == len(self.levels) - 1)))
bclabel = s.get_bclabel()
bstate = s.get_block_state(b=bclabel,
**dict(self.hstate_args,
Lrecdx=s.Lrecdx))
S += bstate.entropy(**dict(self.hentropy_args,
edges_dl=(l + 1 == len(self.levels) - 1)))
if s.Lrecdx[0] >= 0:
S += self._Lrecdx_entropy(s.Lrecdx)
ss = s
while ss is not None:
B_E_D = ss._state.get_B_E_D()
if B_E_D > 0:
for i in range(len(s.rec)):
if s.rec_types[i] != libinference.rec_type.real_normal:
continue
S -= log(B_E_D)
if l < len(self.levels) - 1 and ss is not bstate:
ss = bstate
else:
ss = None
return S
entropy_args = dict(entropy_args, callback=callback)
mcmc_args = dict(mcmc_args, entropy_args=entropy_args)
......
......@@ -243,6 +243,11 @@ class OverlapBlockState(BlockState):
self.Lrecdx = libcore.Vector_double(len(self.rec)+1)
self.Lrecdx[0] = -1
self.Lrecdx.resize(len(self.rec)+1)
self.epsilon = libcore.Vector_double(len(self.rec))
for i in range(len(self.rec)):
idx = self.rec[i].a != 0
if numpy.any(idx):
self.epsilon[i] = abs(self.rec[i].a[idx]).min() / 10
self.max_BE = max_BE
......
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