Commit 253a289f authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

inference.blockmodel: Fix issue with nested SBM + signed edge covariates

parent dc067716
......@@ -356,7 +356,11 @@ public:
}
if (Add && (mrs < 2))
{
if (_B_E_D == 0 && this->_Lrecdx[0] >= 0)
this->_Lrecdx[0] += 1;
_B_E_D++;
}
}
if (mrs > 1)
......@@ -374,7 +378,11 @@ public:
}
if (!Add && ((mrs + get<1>(delta)[0]) < 2))
{
_B_E_D--;
if (_B_E_D == 0 && this->_Lrecdx[0] >= 0)
this->_Lrecdx[0] -= 1;
}
}
for (size_t i = 0; i < this->_rec_types.size(); ++i)
......@@ -440,7 +448,11 @@ public:
_B_E++;
if (_brec[0][me] == 2)
{
if (_B_E_D == 0 && _Lrecdx[0] >= 0)
_Lrecdx[0] += 1;
_B_E_D++;
}
if (_rt == weight_type::REAL_NORMAL && _brec[0][me] > 1)
{
......@@ -482,7 +494,11 @@ public:
_B_E--;
if (_brec[0][me] == 1)
{
_B_E_D--;
if (_B_E_D == 0 && _Lrecdx[0] >= 0)
_Lrecdx[0] -= 1;
}
if (_rt == weight_type::REAL_NORMAL && _brec[0][me] > 0)
{
......@@ -1364,6 +1380,7 @@ public:
}
}
int dL = 0;
if (ea.recs)
{
auto positive_entries_op = [&](size_t i, auto&& w_log_P,
......@@ -1540,6 +1557,24 @@ public:
_recdx[i] + _dBdx[i],
wp[2], wp[3]);
}
if (dL == 0)
{
if (_B_E_D == 0 && dB_E_D > 0)
dL++;
if (_B_E_D > 0 && _B_E_D + dB_E_D == 0)
dL--;
}
if (_coupled_state == nullptr && _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]);
}
}
}
break;
......@@ -1611,10 +1646,11 @@ public:
});
#pragma omp critical (coupled_virtual_move)
{
dS += _coupled_state->recs_dS(r, nr, recs_entries, _dBdx);
dS += _coupled_state->recs_dS(r, nr, recs_entries, _dBdx, dL);
}
}
}
return dS;
}
......@@ -1625,9 +1661,9 @@ public:
double recs_dS(size_t u, size_t v,
const std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
std::vector<double>>>& entries,
std::vector<double>& dBdx)
GraphInterface::edge_t, int,
std::vector<double>>>& entries,
std::vector<double>& dBdx, int dL)
{
if (_rec_types.empty())
return 0;
......@@ -1759,19 +1795,26 @@ public:
std::pow(this->_brec[i][me], 2) / (ers + d));
}
});
if (dB_E_D != 0 || drecdx != 0 || dBdx[i] != 0)
if (dB_E_D != 0 || drecdx != 0 || dBdx[i] != 0 || dL != 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] + drecdx,
wp[2], wp[3]);
size_t L = _Lrecdx[0];
if (L > 0)
int ddL = 0;
if (_Lrecdx[0] >= 0)
{
size_t L = _Lrecdx[0];
if (_B_E_D == 0 && dB_E_D > 0)
ddL++;
if (_B_E_D > 0 && _B_E_D + dB_E_D == 0)
ddL--;
dS -= -positive_w_log_P_alt(L, _Lrecdx[i+1],
wp[2], wp[3]);
dS += -positive_w_log_P_alt(L, _Lrecdx[i+1] +
dS += -positive_w_log_P_alt(L + dL + ddL,
_Lrecdx[i+1] +
dBdx[i] + drecdx,
wp[2], wp[3]);
}
......
......@@ -85,6 +85,10 @@ void export_lsbm()
.def("get_partition_dl", &state_t::get_partition_dl)
.def("get_deg_dl", &state_t::get_deg_dl)
.def("get_move_prob", get_move_prob)
.def("get_B_E",
&state_t::get_B_E)
.def("get_B_E_D",
&state_t::get_B_E_D)
.def("enable_partition_stats",
&state_t::enable_partition_stats)
.def("disable_partition_stats",
......
......@@ -98,6 +98,10 @@ void export_layered_overlap_blockmodel_state()
.def("get_partition_dl", &state_t::get_partition_dl)
.def("get_deg_dl", &state_t::get_deg_dl)
.def("get_move_prob", get_move_prob)
.def("get_B_E",
&state_t::get_B_E)
.def("get_B_E_D",
&state_t::get_B_E_D)
.def("enable_partition_stats",
&state_t::enable_partition_stats)
.def("disable_partition_stats",
......
......@@ -183,6 +183,10 @@ void export_overlap_blockmodel_state()
.def("get_partition_dl", &state_t::get_partition_dl)
.def("get_deg_dl", &state_t::get_deg_dl)
.def("get_move_prob", get_move_prob)
.def("get_B_E",
&state_t::get_B_E)
.def("get_B_E_D",
&state_t::get_B_E_D)
.def("enable_partition_stats",
&state_t::enable_partition_stats)
.def("disable_partition_stats",
......
......@@ -389,6 +389,16 @@ public:
}
}
size_t get_B_E()
{
return _B_E;
}
size_t get_B_E_D()
{
return _B_E_D;
}
void remove_vertex(size_t v)
{
modify_vertex<false>(v, _b[v]);
......@@ -596,6 +606,7 @@ public:
}
}
int dL = 0;
if (ea.recs)
{
auto positive_entries_op = [&](size_t i, auto&& w_log_P,
......@@ -772,6 +783,24 @@ public:
_recdx[i] + _dBdx[i],
wp[2], wp[3]);
}
if (dL == 0)
{
if (_B_E_D == 0 && dB_E_D > 0)
dL++;
if (_B_E_D > 0 && _B_E_D + dB_E_D == 0)
dL--;
}
if (_coupled_state == nullptr && _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]);
}
}
}
break;
......@@ -821,7 +850,8 @@ public:
});
#pragma omp critical (coupled_virtual_move)
{
dS += _coupled_state->recs_dS(r, nr, recs_entries, _dBdx);
dS += _coupled_state->recs_dS(r, nr, recs_entries, _dBdx,
dL);
}
}
}
......@@ -840,9 +870,11 @@ public:
return get_partition_stats(v).get_delta_partition_dl(v, r, nr, _g);
}
double recs_dS(size_t, size_t, const std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int, std::vector<double>>> &,
std::vector<double>&)
double recs_dS(size_t, size_t,
const std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
std::vector<double>>> &,
std::vector<double>&, int)
{
return 0;
}
......
......@@ -1630,7 +1630,7 @@ public:
const std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
std::vector<double>>>&,
std::vector<double>&) = 0;
std::vector<double>&, int) = 0;
};
......
......@@ -429,6 +429,7 @@ class BlockState(object):
self.Lrecdx = kwargs.pop("Lrecdx", None)
if self.Lrecdx is None:
self.Lrecdx = libcore.Vector_double(len(self.rec)+1)
self.Lrecdx[0] = -1
self.Lrecdx.resize(len(self.rec)+1)
self.allow_empty = allow_empty
......@@ -632,7 +633,8 @@ class BlockState(object):
if self._coupled_state is not None:
state._couple_state(state.get_block_state(b=state.get_bclabel(),
copy_bg=False),
copy_bg=False,
Lrecdx=state.Lrecdx),
self._coupled_state[1])
return state
......@@ -775,7 +777,8 @@ class BlockState(object):
if copy_coupled and self._coupled_state is not None:
state._couple_state(state.get_block_state(b=state.get_bclabel(),
copy_bg=False),
copy_bg=False,
Lrecdx=state.Lrecdx),
self._coupled_state[1])
return state
......@@ -1101,7 +1104,8 @@ class BlockState(object):
Salt = state_copy.entropy(**args)
assert math.isclose(S, Salt, abs_tol=1e-8), \
"entropy discrepancy after copying (%g %g)" % (S, Salt)
"entropy discrepancy after copying (%g %g %g)" % (S, Salt,
S - Salt)
if len(kwargs) > 0:
raise ValueError("unrecognized keyword arguments: " +
......
......@@ -35,7 +35,8 @@ import copy
def get_edges_dl(state, hstate_args, hentropy_args):
bclabel = state.get_bclabel()
bstate = state.get_block_state(b=bclabel, **hstate_args)
bstate = state.get_block_state(b=bclabel, **dict(hstate_args,
Lrecdx=state.Lrecdx))
return bstate.entropy(**hentropy_args)
class NestedBlockState(object):
......@@ -111,7 +112,7 @@ class NestedBlockState(object):
def _regen_Lrecdx(self):
self.Lrecdx.a = 0
self.Lrecdx[0] = len(self.levels)
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
......@@ -276,6 +277,15 @@ class NestedBlockState(object):
edges_dl=(l == (len(self.levels) - 1))))
return S
def _Lrecdx_entropy(self, Lrecdx=None):
if Lrecdx is None:
Lrecdx = self.Lrecdx
S = 0
for i in range(len(Lrecdx)-1):
S += -libinference.positive_w_log_P(Lrecdx[0], Lrecdx[i+1],
numpy.nan, numpy.nan)
return S
def entropy(self, **kwargs):
"""Compute the entropy of whole hierarchy.
......@@ -288,10 +298,7 @@ class NestedBlockState(object):
for l in range(len(self.levels)):
S += self.level_entropy(l, **dict(kwargs, test=False))
for i in range(len(self.Lrecdx)-1):
S += -libinference.positive_w_log_P(len(self.levels),
self.Lrecdx[i+1], numpy.nan,
numpy.nan);
S += self._Lrecdx_entropy()
if _bm_test() and kwargs.pop("test", True):
state = self.copy()
......@@ -419,8 +426,8 @@ class NestedBlockState(object):
print("l: %d, N: %d, B: %d" % (l, state.get_N(),
state.get_nonempty_B()))
def find_new_level(self, l, sparse_thres=numpy.inf, bisection_args={}, B_min=None,
B_max=None, b_min=None, b_max=None):
def find_new_level(self, l, sparse_thres=numpy.inf, bisection_args={},
B_min=None, B_max=None, b_min=None, b_max=None):
"""Attempt to find a better network partition at level ``l``, using
:func:`~graph_tool.inference.bisection_minimize` with arguments given by
``bisection_args``.
......@@ -444,13 +451,16 @@ class NestedBlockState(object):
else:
entropy_args = dict(entropy_args, dl=True,
edges_dl=l==len(self.levels) - 1)
if l < len(self.levels) - 1:
entropy_args = dict(entropy_args,
callback=\
lambda s: get_edges_dl(s,
self.hstate_args,
dict(self.hentropy_args,
edges_dl=(l + 1 == len(self.levels) - 1))))
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)))
if s.Lrecdx[0] >= 0:
S += self._Lrecdx_entropy(s.Lrecdx)
return S
entropy_args = dict(entropy_args, callback=callback)
mcmc_args = dict(mcmc_args, entropy_args=entropy_args)
if _bm_test() and isinstance(self.levels[0], LayeredBlockState):
mcmc_args = dict(mcmc_args, disable_callback_test=True)
......@@ -871,6 +881,7 @@ def hierarchy_minimize(state, B_min=None, B_max=None, b_min=None, b_max=None,
Sf - Si, len(state.levels))
else:
state.levels[l:l+len(bstates)] = bstates
state._regen_Lrecdx()
if check_verbose(verbose):
print(verbose_pad(verbose) + "level", l,
......@@ -894,6 +905,7 @@ def hierarchy_minimize(state, B_min=None, B_max=None, b_min=None, b_max=None,
if Si - Sf < epsilon:
state.levels[l - 1] = bstates[0]
state.levels.insert(l, bstates[1])
state._regen_Lrecdx()
else:
kept = False
del done[l]
......@@ -936,6 +948,7 @@ def hierarchy_minimize(state, B_min=None, B_max=None, b_min=None, b_max=None,
state.levels[l + j] = bstates[j]
if bstates[-1].B == 1:
del state.levels[l + len(bstates):]
state._regen_Lrecdx()
else:
kept = False
dS += Sf - Si
......@@ -953,6 +966,8 @@ def hierarchy_minimize(state, B_min=None, B_max=None, b_min=None, b_max=None,
bstate = bstate.get_block_state(b=zeros(state.levels[-1].B),
deg_corr=False)
state.levels.append(bstate)
state._regen_Lrecdx()
if _bm_test():
state._consistency_check()
......
......@@ -238,8 +238,10 @@ class OverlapBlockState(BlockState):
BlockState._init_recs(self, self.rec, rec_types, rec_params)
self.recdx = libcore.Vector_double(len(self.rec))
self.Lrecdx = kwargs.pop("Lrecdx",
libcore.Vector_double(len(self.rec)+1))
self.Lrecdx = kwargs.pop("Lrecdx", None)
if self.Lrecdx is None:
self.Lrecdx = libcore.Vector_double(len(self.rec)+1)
self.Lrecdx[0] = -1
self.Lrecdx.resize(len(self.rec)+1)
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