Commit 878fb5b6 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

inference: remove 'allow_empty' option

parent 67372fc8
Pipeline #520 failed with stage
in 622 minutes and 14 seconds
......@@ -38,7 +38,7 @@ rec_p = g.new_ep("double", random(g.num_edges()))
rec_s = g.new_ep("double", normal(0, 10, g.num_edges()))
def _gen_state(directed, deg_corr, layers, overlap, rec, rec_type, allow_empty):
def _gen_state(directed, deg_corr, layers, overlap, rec, rec_type):
u = GraphView(g, directed=directed)
if layers != False:
base_type = graph_tool.inference.LayeredBlockState
......@@ -48,8 +48,7 @@ def _gen_state(directed, deg_corr, layers, overlap, rec, rec_type, allow_empty):
recs=[rec] if rec is not None else [],
rec_types=[rec_type] if rec is not None else [],
overlap=overlap,
layers=layers == True,
allow_empty=allow_empty)
layers=layers == True)
elif overlap:
base_type = graph_tool.inference.OverlapBlockState
state_args = dict(B=2 * u.num_edges(),
......@@ -61,8 +60,7 @@ def _gen_state(directed, deg_corr, layers, overlap, rec, rec_type, allow_empty):
state_args = dict(B=u.num_vertices(),
recs=[rec] if rec is not None else [],
rec_types=[rec_type] if rec is not None else [],
deg_corr=deg_corr,
allow_empty=allow_empty)
deg_corr=deg_corr)
return u, base_type, state_args
......@@ -71,9 +69,8 @@ def gen_state(*args):
return base_type(u, **state_args)
def gen_nested_state(*args):
u, base_type, state_args = _gen_state(*args, False)
u, base_type, state_args = _gen_state(*args)
B = state_args.pop("B")
state_args.pop("allow_empty", None)
return NestedBlockState(u,
bs=[numpy.arange(B)] + [numpy.zeros(1)] * 6,
base_type=base_type, state_args=state_args,
......@@ -87,7 +84,6 @@ pranges = [("directed", [False, True]),
("deg_corr", [False, True]),
("dl", [False, True]),
("degree_dl_kind", ["uniform", "distributed", "entropy"]),
("allow_empty", [False, True]),
("exact", [True, False])]
pranges = OrderedDict(pranges)
......@@ -120,7 +116,7 @@ for pvals in iter_ranges(pranges):
rec_ = rec_s
print("\t mcmc (unweighted)", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
print("\t\t",
state.mcmc_sweep(beta=0,
......@@ -138,7 +134,7 @@ for pvals in iter_ranges(pranges):
state.get_nonempty_B(), file=out)
print("\t mcmc (unweighted, multiflip)", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
print("\t\t",
state.multiflip_mcmc_sweep(beta=0,
entropy_args=dict(dl=dl,
......@@ -147,7 +143,7 @@ for pvals in iter_ranges(pranges):
state.get_nonempty_B(), file=out)
print("\t gibbs (unweighted)", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
print("\t\t",
state.gibbs_sweep(beta=0,
......@@ -157,7 +153,7 @@ for pvals in iter_ranges(pranges):
state.get_nonempty_B(), file=out)
if not overlap:
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
print("\t mcmc", file=out)
bstate = state.get_block_state(vweight=True, deg_corr=deg_corr)
......@@ -185,7 +181,7 @@ for pvals in iter_ranges(pranges):
print("\t merge", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
if not overlap:
bstate = state.get_block_state(vweight=True, deg_corr=deg_corr)
......@@ -223,7 +219,7 @@ for pvals in iter_ranges(pranges):
print("\t shrink", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec)
state = state.shrink(B=5, entropy_args=dict(dl=dl,
degree_dl_kind=degree_dl_kind,
multigraph=False,
......
......@@ -101,8 +101,7 @@ typedef mpl::vector1<std::false_type> rmap_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))
((epsilon, &, std::vector<double>&, 0))
GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)
......@@ -271,7 +270,7 @@ public:
BlockState::remove_partition_node(v, r);
}
bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true)
bool allow_move(size_t v, size_t r, size_t nr)
{
if (_coupled_state != nullptr && is_last(v))
{
......@@ -280,10 +279,7 @@ public:
return false;
}
if (allow_empty)
return ((_bclabel[r] == _bclabel[nr]) || (_wr[nr] == 0));
else
return _bclabel[r] == _bclabel[nr];
return _bclabel[r] == _bclabel[nr];
}
template <class EFilt>
......@@ -1665,8 +1661,6 @@ public:
enable_partition_stats();
for (auto& ps : _partition_stats)
actual_B += ps.get_actual_B();
if (_allow_empty)
actual_B = num_vertices(_bg);
size_t E = _partition_stats.front().get_E();
dS -= get_edges_dl(actual_B, E, _g);
dS += get_edges_dl(actual_B + dr + ds, E, _g);
......@@ -2051,8 +2045,6 @@ public:
size_t actual_B = 0;
for (auto& ps : _partition_stats)
actual_B += ps.get_actual_B();
if (_allow_empty)
actual_B = num_vertices(_bg);
S_dl += get_edges_dl(actual_B, _partition_stats.front().get_E(), _g);
}
......@@ -2381,7 +2373,7 @@ public:
for (size_t c = 0; c < C; ++c)
_partition_stats.emplace_back(_g, _b, vcs[c], E, B,
_vweight, _eweight, _degs,
_bmap, _allow_empty);
_bmap);
for (auto r : vertices_range(_bg))
_partition_stats[rc[r]].get_r(r);
......
......@@ -95,7 +95,7 @@ struct Gibbs
{
if (nr == null_group)
{
if (!_allow_new_group || _state._allow_empty)
if (!_allow_new_group)
return numeric_limits<double>::infinity();
if (_state._empty_blocks.empty())
_state.add_block();
......
......@@ -113,7 +113,7 @@ struct Merge
s = uniform_sample(_available, rng);
}
if (s == r || !_state.allow_move(v, r, s, false))
if (s == r || !_state.allow_move(v, r, s))
return _null_move;
return s;
......
......@@ -71,8 +71,8 @@ public:
class Vlist>
partition_stats(Graph& g, Vprop& b, Vlist& vlist, size_t E, size_t B,
VWprop& vweight, Eprop& eweight, Degs& degs,
std::vector<size_t>& bmap, bool allow_empty)
: _bmap(bmap), _N(0), _E(E), _total_B(B), _allow_empty(allow_empty)
std::vector<size_t>& bmap)
: _bmap(bmap), _N(0), _E(E), _total_B(B)
{
if (!use_rmap)
{
......@@ -134,10 +134,7 @@ public:
double get_partition_dl()
{
double S = 0;
if (_allow_empty)
S += lbinom(_total_B + _N - 1, _N);
else
S += lbinom(_N - 1, _actual_B - 1);
S += lbinom(_N - 1, _actual_B - 1);
S += lgamma_fast(_N + 1);
for (auto nr : _total)
S -= lgamma_fast(nr + 1);
......@@ -304,7 +301,7 @@ public:
if (nr != null_group && _total[nr] == 0)
dB++;
if ((dN != 0 || dB != 0) && !_allow_empty)
if ((dN != 0 || dB != 0))
{
S_b += lbinom_fast(_N - 1, _actual_B - 1);
S_a += lbinom_fast(_N - 1 + dN, _actual_B + dB - 1);
......@@ -323,7 +320,7 @@ public:
double get_delta_edges_dl(size_t v, size_t r, size_t nr, VProp& vweight,
size_t actual_B, Graph& g)
{
if (r == nr || _allow_empty)
if (r == nr)
return 0;
if (r != null_group)
......@@ -670,7 +667,6 @@ private:
size_t _E;
size_t _actual_B;
size_t _total_B;
bool _allow_empty;
vector<map_t> _hist;
vector<int> _total;
vector<int> _ep;
......
......@@ -71,7 +71,7 @@ public:
entropy_args_t& ea) = 0;
virtual vprop_map_t<int32_t>::type::unchecked_t& get_b() = 0;
virtual bool check_edge_counts(bool emat=true) = 0;
virtual bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true) = 0;
virtual bool allow_move(size_t v, size_t r, size_t nr) = 0;
};
} // graph_tool namespace
......
......@@ -95,13 +95,11 @@ struct Layers
{
public:
LayerState(const BaseState& base_state, LayeredBlockState& lstate,
bmap_t& block_map, block_rmap_t block_rmap,
vector<size_t>& free_blocks, size_t l)
bmap_t& block_map, block_rmap_t block_rmap, size_t l)
: BaseState(base_state),
_lstate(&lstate),
_block_map(block_map),
_block_rmap(block_rmap),
_free_blocks(free_blocks),
_l(l), _E(0)
{
for (auto e : edges_range(BaseState::_g))
......@@ -111,12 +109,12 @@ struct Layers
LayeredBlockState* _lstate;
bmap_t& _block_map;
block_rmap_t _block_rmap;
vector<size_t>& _free_blocks;
size_t _l;
size_t _E;
using BaseState::_bg;
using BaseState::_wr;
using BaseState::_empty_blocks;
using BaseState::add_block;
size_t get_block_map(size_t r, bool put_new=true)
......@@ -125,11 +123,19 @@ struct Layers
auto iter = _block_map.find(r);
if (iter == _block_map.end())
{
if (_free_blocks.empty())
_free_blocks.push_back(_block_map.size());
r_u = _free_blocks.back();
while (r_u >= num_vertices(_bg))
add_block();
r_u = null_group;
for (auto s : _empty_blocks)
{
if (_block_rmap[s] != -1)
continue;
r_u = s;
break;
}
if (r_u == null_group)
{
r_u = add_block();
_block_rmap[r_u] = -1;
}
assert(r_u < num_vertices(_bg));
if (put_new)
......@@ -147,7 +153,6 @@ struct Layers
// _lstate->_lcoupled_state->set_block(_l, r_u, hb[r_u]);
// assert(_lstate->_lcoupled_state->get_block(_l, r_u) == size_t(hb[r_u]));
}
_free_blocks.pop_back();
assert(_lstate->_lcoupled_state == nullptr ||
r_u == _lstate->_lcoupled_state->get_layer_node(_l, r));
// assert(_lstate->_lcoupled_state == nullptr ||
......@@ -204,10 +209,8 @@ struct Layers
boost::python::object temp = ostate.attr("block_rmap").attr("_get_any")();
boost::any& a = python::extract<boost::any&>(temp);
block_rmap_t block_rmap = boost::any_cast<block_rmap_t>(a);
std::vector<size_t>& free_blocks =
python::extract<std::vector<size_t>&>(ostate.attr("free_blocks"));
bmap_t& block_map = _block_map[l];
_layers.emplace_back(state, *this, block_map, block_rmap, free_blocks, l);
_layers.emplace_back(state, *this, block_map, block_rmap, l);
}
for (auto r : vertices_range(BaseState::_bg))
{
......@@ -487,9 +490,9 @@ struct Layers
set_partition(b.get_unchecked());
}
bool allow_move(size_t v, size_t r, size_t nr, bool allow_empty = true)
bool allow_move(size_t v, size_t r, size_t nr)
{
return BaseState::allow_move(v, r, nr, allow_empty);
return BaseState::allow_move(v, r, nr);
}
template <class MEntries>
......@@ -695,11 +698,8 @@ struct Layers
if (ea.edges_dl)
{
size_t actual_B = _actual_B;
if (BaseState::_allow_empty)
actual_B = num_vertices(BaseState::_bg);
for (auto& state : _layers)
S_dl += get_edges_dl(actual_B, state._E, _g);
S_dl += get_edges_dl(_actual_B, state._E, _g);
}
if (ea.recs)
......@@ -729,16 +729,9 @@ struct Layers
for (auto& state : _layers)
{
size_t actual_B = 0;
if (BaseState::_allow_empty)
{
actual_B = num_vertices(_bg);
}
else
{
for (auto r : vertices_range(state._bg))
if (state._wr[r] > 0)
actual_B++;
}
for (auto r : vertices_range(state._bg))
if (state._wr[r] > 0)
actual_B++;
S_dl += get_edges_dl(actual_B, state._E, _g);
}
}
......@@ -754,7 +747,7 @@ struct Layers
double get_delta_edges_dl(size_t v, size_t r, size_t s)
{
if (r == s || BaseState::_allow_empty)
if (r == s)
return 0;
if (BaseState::_vweight[v] == 0)
return 0;
......@@ -768,7 +761,7 @@ struct Layers
{
auto get_x = [](size_t B)
{
if (is_directed_::apply<typename BaseState::g_t>::type::value)
if constexpr (is_directed_::apply<typename BaseState::g_t>::type::value)
return B * B;
else
return (B * (B + 1)) / 2;
......@@ -957,15 +950,15 @@ struct Layers
size_t add_block()
{
auto r = BaseState::add_block();
for (size_t l = 0; l < _layers.size(); ++l)
{
auto& state = _layers[l];
size_t r_u = state.add_block();
if (_lcoupled_state != nullptr)
_lcoupled_state->get_layer(l).coupled_resize_vertex(r_u);
}
return r;
return BaseState::add_block();
// for (size_t l = 0; l < _layers.size(); ++l)
// {
// auto& state = _layers[l];
// size_t r_u = state.add_block();
// if (_lcoupled_state != nullptr)
// _lcoupled_state->get_layer(l).coupled_resize_vertex(r_u);
// }
// return r;
}
void coupled_resize_vertex(size_t v)
......
......@@ -50,7 +50,10 @@ typedef mpl::vector2<std::true_type, std::false_type> use_hash_tr;
((mrm,, vmap_t, 0)) \
((wr,, vmap_t, 0)) \
((b,, vmap_t, 0)) \
((empty_blocks, & ,std::vector<size_t>&, 0)) \
((empty_pos,, vmap_t, 0)) \
((candidate_blocks, &, std::vector<size_t>&, 0)) \
((candidate_pos,, vmap_t, 0)) \
((bclabel,, vmap_t, 0)) \
((pclabel,, vmap_t, 0)) \
((deg_corr,, bool, 0)) \
......@@ -63,8 +66,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))
((epsilon, &, std::vector<double>&, 0))
GEN_STATE_BASE(OverlapBlockStateVirtualBase, OVERLAP_BLOCK_STATE_params)
......@@ -88,8 +90,17 @@ public:
_overlap_stats(_g, _b, _half_edges, _node_index, num_vertices(_bg)),
_coupled_state(nullptr)
{
_empty_blocks.clear();
_candidate_blocks.clear();
for (auto r : vertices_range(_bg))
{
_wr[r] = _overlap_stats.get_block_size(r);
if (_wr[r] == 0)
add_element(_empty_blocks, _empty_pos, r);
else
add_element(_candidate_blocks, _candidate_pos, r);
}
for (auto& p : _brec)
{
_c_brec.push_back(p.get_checked());
......@@ -163,6 +174,12 @@ public:
template <bool Add>
void modify_vertex(size_t v, size_t r)
{
if (Add && _wr[r] == 0)
{
remove_element(_empty_blocks, _empty_pos, r);
add_element(_candidate_blocks, _candidate_pos, r);
}
if (Add)
get_move_entries(v, null_group, r, _m_entries);
else
......@@ -183,7 +200,14 @@ public:
if (!_egroups.empty() && _egroups_enabled)
_egroups.remove_vertex(v, _b, _g);
}
_wr[r] = _overlap_stats.get_block_size(r);
if (!Add && _wr[r] == 0)
{
remove_element(_candidate_blocks, _candidate_pos, r);
add_element(_empty_blocks, _empty_pos, r);
}
}
size_t get_B_E()
......@@ -206,14 +230,16 @@ public:
modify_vertex<true>(v, r);
}
bool allow_move(size_t, size_t r, size_t nr)
bool allow_move(size_t v, size_t r, size_t nr)
{
return (_bclabel[r] == _bclabel[nr]);
}
if (_coupled_state != nullptr && is_last(v))
{
auto& bh = _coupled_state->get_b();
if (bh[r] != bh[nr])
return false;
}
bool allow_move(size_t v, size_t r, size_t nr, bool)
{
return allow_move(v, r, nr);
return _bclabel[r] == _bclabel[nr];
}
// move a vertex from its current block to block nr
......@@ -508,10 +534,10 @@ public:
{
if (r_vacate)
dS += _coupled_state->get_delta_partition_dl(r, bh[r], null_group,
_coupled_entropy_args);
_coupled_entropy_args);
if (nr_occupy)
dS += _coupled_state->get_delta_partition_dl(nr, null_group, bh[nr],
_coupled_entropy_args);
_coupled_entropy_args);
}
}
return dS;
......@@ -519,11 +545,28 @@ public:
// Sample node placement
template <class RNG>
size_t sample_block(size_t v, double c, double, RNG& rng)
size_t sample_block(size_t v, double c, double d, RNG& rng)
{
// attempt random block
size_t s = uniform_sample(_candidate_blocks, rng);
// attempt new block
std::bernoulli_distribution new_r(d);
if (d > 0 && new_r(rng) && (_candidate_blocks.size() < num_vertices(_g)))
{
if (_empty_blocks.empty())
add_block();
s = uniform_sample(_empty_blocks, rng);
auto r = _b[v];
_bclabel[s] = _bclabel[r];
if (_coupled_state != nullptr)
{
auto& hb = _coupled_state->get_b();
hb[s] = hb[r];
}
return s;
}
if (!std::isinf(c))
{
size_t w = get_lateral_half_edge(v, rng);
......@@ -536,7 +579,7 @@ public:
double p_rand = 0;
if (c > 0)
{
size_t B = num_vertices(_bg);
size_t B = _candidate_blocks.size();
if (graph_tool::is_directed(_g))
p_rand = c * B / double(_mrp[t] + _mrm[t] + c * B);
else
......@@ -555,16 +598,6 @@ public:
}
}
// if (_wr[s] == 0 && _coupled_state != nullptr)
// {
// auto& hb = _coupled_state->get_b();
// do
// {
// hb[s] = _coupled_state->sample_block(_b[v], std::numeric_limits<double>::infinity(), d, rng);
// }
// while (!_coupled_state->allow_move(hb[_b[v]], hb[s], _allow_empty));
// }
return s;
}
......@@ -594,11 +627,32 @@ public:
// Computes the move proposal probability
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c, double,
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse, MEntries& m_entries)
{
size_t B = _candidate_blocks.size();
if (reverse)
{
if (_overlap_stats.virtual_remove_size(v, s) == 0)
return d;
if (_wr[r] == 0)
B++;
}
else
{
if (_wr[s] == 0)
return d;
}
if (B == num_vertices(_g))
d = 0;
if (std::isinf(c))
return (1. - d) / B;
typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
size_t B = num_vertices(_bg);
double p = 0;
size_t w = 0;
......@@ -676,9 +730,9 @@ public:
}
}
if (w > 0)
return p / w;
return (1. - d) * p / w;
else
return 1. / B;
return (1. - d) / B;
}
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
......@@ -793,8 +847,6 @@ public:
size_t actual_B = 0;
for (auto& ps : _partition_stats)
actual_B += ps.get_actual_B();