Skip to content
GitLab
Commit 4f0c7fe6 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

inference.blockmodel: Improve nested SBM equilibration 

This modifies the MCMC code to allow simultaneous update of
all hierarchy levels, improving equilibration speed.
parent 749572d8
Expand all Hide whitespace changes
Inline Side-by-side
doc/test_inference.py
View file @ 4f0c7fe6
......@@ -23,12 +23,12 @@ graph_tool.inference.set_test(True)
g = collection.data["football"]
# add self-loops
# # add self-loops
# for i in range(10):
# v = numpy.random.randint(g.num_vertices())
# g.add_edge(v, v)
# add parallel edges
# # add parallel edges
# for e in list(g.edges())[:10]:
# g.add_edge(e.source(), e.target())
......@@ -38,29 +38,46 @@ 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, allow_empty):
u = GraphView(g, directed=directed)
if layers != False:
state = graph_tool.inference.LayeredBlockState(u, B=u.num_vertices(),
deg_corr=deg_corr,
ec=ec.copy(),
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)
base_type = graph_tool.inference.LayeredBlockState
state_args = dict(B=u.num_vertices(),
deg_corr=deg_corr,
ec=ec.copy(),
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)
elif overlap:
state = graph_tool.inference.OverlapBlockState(u, B=2 * u.num_edges(),
recs=[rec] if rec is not None else [],
rec_types=[rec_type] if rec is not None else [],
deg_corr=deg_corr)
base_type = graph_tool.inference.OverlapBlockState
state_args = dict(B=2 * u.num_edges(),
recs=[rec] if rec is not None else [],
rec_types=[rec_type] if rec is not None else [],
deg_corr=deg_corr)
else:
state = graph_tool.inference.BlockState(u, 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)
return state
base_type = graph_tool.inference.BlockState
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)
return u, base_type, state_args
def gen_state(*args):
u, base_type, state_args = _gen_state(*args)
return base_type(u, **state_args)
def gen_nested_state(*args):
u, base_type, state_args = _gen_state(*args, False)
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,
sampling=True)
pranges = [("directed", [False, True]),
......@@ -120,9 +137,28 @@ for pvals in iter_ranges(pranges):
exact=exact, beta_dl=0.95)),
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)
print("\t\t",
state.multiflip_mcmc_sweep(beta=0,
entropy_args=dict(dl=dl,
degree_dl_kind=degree_dl_kind,
exact=exact, beta_dl=0.95)),
state.get_nonempty_B(), file=out)
print("\t gibbs (unweighted)", file=out)
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
print("\t\t",
state.gibbs_sweep(beta=0,
entropy_args=dict(dl=dl,
degree_dl_kind=degree_dl_kind,
exact=exact, beta_dl=0.95)),
state.get_nonempty_B(), file=out)
if not overlap:
state = gen_state(directed, deg_corr, layered, overlap, rec_, rec, allow_empty)
print("\t mcmc", file=out)
bstate = state.get_block_state(vweight=True, deg_corr=deg_corr)
......@@ -194,6 +230,59 @@ for pvals in iter_ranges(pranges):
exact=exact, beta_dl=0.95))
print("\t\t", state.B, "\n", file=out)
pranges = [("directed", [False, True]),
("overlap", [False]),
("layered", [False, "covariates", True]),
("rec", [None, "real-exponential", "real-normal"]),
("deg_corr", [True, False]),
("degree_dl_kind", ["distributed"]),
("exact", [True])]
pranges = OrderedDict(pranges)
def iter_ranges(ranges):
for vals in itertools.product(*[v for k, v in ranges.items()]):
yield zip(ranges.keys(), vals)
for pvals in iter_ranges(pranges):
params = OrderedDict(pvals)
locals().update(params)
if overlap and deg_corr and degree_dl_kind != "distributed": # FIXME
continue
if (rec is not None or layered != False) and not exact:
continue
print(params, file=out)
rec_ = None
if rec == "real-exponential":
rec_ = rec_p
elif rec == "real-normal":
rec_ = rec_s
print("\t mcmc (single flip)", file=out)
state = gen_nested_state(directed, deg_corr, layered, overlap, rec_, rec)
for i in range(5):
print("\t\t",
state.mcmc_sweep(beta=0, d=0.5,
entropy_args=dict(degree_dl_kind=degree_dl_kind,
exact=exact,
beta_dl=0.95)),
[s.get_nonempty_B() for s in state.levels], file=out)
print("\n\t mcmc (multiple flip)", file=out)
state = gen_nested_state(directed, deg_corr, layered, overlap, rec_, rec)
for i in range(5):
print("\t\t",
state.multiflip_mcmc_sweep(beta=0, d=0.5,
entropy_args=dict(degree_dl_kind=degree_dl_kind,
exact=exact, beta_dl=0.95)),
[s.get_nonempty_B() for s in state.levels], file=out)
pranges = [("directed", [False, True]),
("overlap", [False, True]),
......
doc/test_inference_mcmc.py
View file @ 4f0c7fe6
......@@ -104,89 +104,107 @@ for i in range(3 * 8):
g.add_edge(s, t)
for directed in [True, False]:
g.set_directed(directed)
hists = {}
for nested in [False, True]:
for directed in [False, True]:
g.set_directed(directed)
state = minimize_blockmodel_dl(g, deg_corr=True)
state = state.copy(B=g.num_vertices())
hists = {}
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs", -1]))
if nested:
state = minimize_nested_blockmodel_dl(g, deg_corr=True)
state = state.copy(bs=[g.get_vertices()] + [zeros(1)] * 6, sampling=True)
else:
state = minimize_blockmodel_dl(g, deg_corr=True)
state = state.copy(B=g.num_vertices())
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=1, c=abs(c), niter=40)
if nested:
cs = list(reversed([-numpy.inf, numpy.inf, 1, 0.1, 0.01, -1]))
else:
mcmc_args=dict(beta=1, niter=40)
if i == 0:
mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
nbreaks=5,
wait=1000,
verbose=(1, "c = %s (t) " % str(c)) if verbose else False)
hists[c] = mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
wait=1000,
nbreaks=5,
verbose=(1, "c = %s " % str(c)) if verbose else False,
history=True)
for c1 in cs:
for c2 in cs:
try:
if c2 < c1:
continue
except TypeError:
pass
Ss1 = array(list(zip(*hists[c1]))[2])
Ss2 = array(list(zip(*hists[c2]))[2])
# add very small normal noise, to solve discreteness issue
Ss1 += numpy.random.normal(0, 1e-2, len(Ss1))
Ss2 += numpy.random.normal(0, 1e-2, len(Ss2))
D, p = scipy.stats.ks_2samp(Ss1, Ss2)
D_c = 1.63 * sqrt((len(Ss1) + len(Ss2)) / (len(Ss1) * len(Ss2)))
if verbose:
print("directed:", directed, "c1:", c1, "c2:", c2,
"D:", D, "D_c:", D_c, "p-value:", p)
if p < .001:
print(("Warning, distributions for directed=%s (c1, c2) = " +
"(%s, %s) are not the same, with a p-value: %g (D=%g, D_c=%g)") %
(str(directed), str(c1), str(c2), p, D, D_c))
for cum in [True, False]:
clf()
bins = None
for c in cs:
hist = hists[c]
if cum:
h = histogram(list(zip(*hist))[0], 1000000, density=True)
y = numpy.cumsum(h[0])
y /= y[-1]
plot(h[-1][:-1], y, "-", label="c=%s" % str(c))
if c != numpy.inf:
hist = hists[numpy.inf]
h2 = histogram(list(zip(*hist))[0], bins=h[-1], density=True)
y2 = numpy.cumsum(h2[0])
y2 /= y2[-1]
res = abs(y - y2)
i = res.argmax()
axvline(h[-1][i], color="grey")
cs = list(reversed([-numpy.inf, numpy.inf, 1, 0.1, 0.01, "gibbs", -1]))
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=.2, c=abs(c), niter=40)
else:
h = histogram(list(zip(*hist))[0], 40, density=True)
plot(h[-1][:-1], h[0], "s-", label="c=%s" % str(c))
gca().set_yscale("log")
mcmc_args=dict(beta=.2, niter=40)
if cum:
legend(loc="lower right")
else:
legend(loc="best")
savefig("test_mcmc_directed%s-cum%s.pdf" % (directed, str(cum)))
if nested:
get_B = lambda s: [sl.get_nonempty_B() for sl in s.levels]
else:
get_B = lambda s: [s.get_nonempty_B()]
if i == 0:
mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
nbreaks=5,
wait=1000,
callback=get_B,
verbose=(1, "c = %s (t) " % str(c)) if verbose else False)
hists[c] = mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
wait=1000,
nbreaks=5,
callback=get_B,
verbose=(1, "c = %s " % str(c)) if verbose else False,
history=True)
output = open("test_mcmc_nested%s_directed%s-output" % (nested, directed), "w")
for c1 in cs:
for c2 in cs:
try:
if c2 < c1:
continue
except TypeError:
pass
Ss1 = array(list(zip(*hists[c1]))[2])
Ss2 = array(list(zip(*hists[c2]))[2])
# add very small normal noise, to solve discreteness issue
Ss1 += numpy.random.normal(0, 1e-2, len(Ss1))
Ss2 += numpy.random.normal(0, 1e-2, len(Ss2))
D, p = scipy.stats.ks_2samp(Ss1, Ss2)
D_c = 1.63 * sqrt((len(Ss1) + len(Ss2)) / (len(Ss1) * len(Ss2)))
print("nested:", nested, "directed:", directed, "c1:", c1,
"c2:", c2, "D:", D, "D_c:", D_c, "p-value:", p,
file=output)
if p < .001:
print(("Warning, distributions for directed=%s (c1, c2) = " +
"(%s, %s) are not the same, with a p-value: %g (D=%g, D_c=%g)") %
(str(directed), str(c1), str(c2), p, D, D_c),
file=output)
output.flush()
for cum in [True, False]:
clf()
bins = None
for c in cs:
hist = hists[c]
if cum:
h = histogram(list(zip(*hist))[2], 1000000, density=True)
y = numpy.cumsum(h[0])
y /= y[-1]
plot(h[-1][:-1], y, "-", label="c=%s" % str(c))
if c != numpy.inf:
hist = hists[numpy.inf]
h2 = histogram(list(zip(*hist))[2], bins=h[-1], density=True)
y2 = numpy.cumsum(h2[0])
y2 /= y2[-1]
res = abs(y - y2)
i = res.argmax()
axvline(h[-1][i], color="grey")
else:
h = histogram(list(zip(*hist))[2], 40, density=True)
plot(h[-1][:-1], h[0], "s-", label="c=%s" % str(c))
gca().set_yscale("log")
if cum:
legend(loc="lower right")
else:
legend(loc="best")
savefig("test_mcmc_nested%s_directed%s-cum%s.pdf" % (nested, directed, str(cum)))
print("OK")
\ No newline at end of file
src/graph/generation/dynamic_sampler.hh
View file @ 4f0c7fe6
......@@ -157,7 +157,10 @@ public:
{
if (_idx[i] == numeric_limits<size_t>::max())
continue;
items.push_back(_items[_idx[i]]);
size_t j = _idx[i];
if (!_valid[j])
continue;
items.push_back(_items[j]);
probs.push_back(_tree[i]);
}
......
src/graph/inference/blockmodel/graph_blockmodel_elist.hh
View file @ 4f0c7fe6
......@@ -133,6 +133,7 @@ public:
template <class Edge, class Vprop>
void insert_edge(const Edge& e, size_t weight, Vprop& b, Graph& g)
{
assert(e != Edge());
size_t r = b[get_source(e, g)];
auto& r_elist = _egroups[r];
insert_edge(std::make_tuple(e, true), r_elist, weight,
......@@ -171,10 +172,10 @@ public:
auto& pos = _epos[e];
size_t r = b[get_source(e, g)];
remove_edge(make_tuple(e, true), pos.first, _egroups[r]);
remove_edge(std::make_tuple(e, true), pos.first, _egroups[r]);
size_t s = b[get_target(e, g)];
remove_edge(make_tuple(e, false), pos.second, _egroups[s]);
remove_edge(std::make_tuple(e, false), pos.second, _egroups[s]);
}
template <class Edge>
......@@ -209,7 +210,7 @@ public:
elist.remove(pos);
if (elist.empty())
elist.rebuild();
elist.clear(true);
}
template <class Vertex, class VProp>
......
src/graph/inference/blockmodel/graph_blockmodel_emat.hh
View file @ 4f0c7fe6
......@@ -55,13 +55,14 @@ public:
typedef typename graph_traits<BGraph>::vertex_descriptor vertex_t;
typedef typename graph_traits<BGraph>::edge_descriptor edge_t;
const auto& get_me(vertex_t r, vertex_t s) const
const edge_t& get_me(vertex_t r, vertex_t s) const
{
return _mat[r][s];
}
void put_me(vertex_t r, vertex_t s, const edge_t& e)
{
assert(e != _null_edge);
_mat[r][s] = e;
if (!is_directed_::apply<BGraph>::type::value && r != s)
_mat[s][r] = e;
......@@ -77,7 +78,7 @@ public:
//remove_edge(me, bg);
}
const auto& get_null_edge() const { return _null_edge; }
const edge_t& get_null_edge() const { return _null_edge; }
private:
multi_array<edge_t, 2> _mat;
......@@ -164,7 +165,7 @@ public:
typedef typename graph_traits<BGraph>::edge_descriptor edge_t;
__attribute__((flatten)) __attribute__((hot))
const auto& get_me(vertex_t r, vertex_t s) const
const edge_t& get_me(vertex_t r, vertex_t s) const
{
if (!is_directed_::apply<BGraph>::type::value && r > s)
std::swap(r, s);
......@@ -180,6 +181,7 @@ public:
if (!is_directed_::apply<BGraph>::type::value && r > s)
std::swap(r, s);
assert(r < _hash.size());
assert(e != _null_edge);
_hash[r][s] = e;
}
......@@ -197,7 +199,7 @@ public:
//remove_edge(me, bg);
}
const auto& get_null_edge() const { return _null_edge; }
const edge_t& get_null_edge() const { return _null_edge; }
private:
std::vector<size_t> _index;
......
src/graph/inference/blockmodel/graph_blockmodel_entries.hh
View file @ 4f0c7fe6
......@@ -18,6 +18,7 @@
#ifndef GRAPH_BLOCKMODEL_ENTRIES_HH
#define GRAPH_BLOCKMODEL_ENTRIES_HH
#include "graph_blockmodel_weights.hh"
namespace std
{
......@@ -246,7 +247,7 @@ public:
_delta.clear();
_edelta.clear();
_mes.clear();
_recs_entries.clear();
_p_entries.clear();
}
const vector<pair<size_t, size_t>>& get_entries() { return _entries; }
......@@ -279,7 +280,7 @@ public:
std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int, std::vector<double>>>
_recs_entries;
_p_entries;
private:
static constexpr size_t _null = numeric_limits<size_t>::max();
......@@ -503,6 +504,90 @@ double entries_dS(MEntries& m_entries, Eprop& mrs, EMat& emat, BGraph& bg)
return dS;
}
template <bool Add, bool Remove, class State, class MEntries>
void apply_delta(State& state, MEntries& m_entries)
{
auto eops =
[&](auto&& eop, auto&& mid_op, auto&& end_op, auto&& skip)
{
eop(m_entries, state._emat,
[&](auto r, auto s, auto& me, auto delta, auto&... edelta)
{
if (skip(delta, edelta...))
return;
if (Add && me == state._emat.get_null_edge())
{
me = boost::add_edge(r, s, state._bg).first;
state._emat.put_me(r, s, me);
state._c_mrs[me] = 0;
for (size_t i = 0; i < state._rec_types.size(); ++i)
{
state._c_brec[i][me] = 0;
state._c_bdrec[i][me] = 0;
}
if (state._coupled_state != nullptr)
state._coupled_state->add_edge(me);
}
mid_op(me, edelta...);
state._mrs[me] += delta;
state._mrp[r] += delta;
state._mrm[s] += delta;
assert(state._mrs[me] >= 0);
assert(state._mrp[r] >= 0);
assert(state._mrm[s] >= 0);
end_op(me, edelta...);
if (Remove && state._mrs[me] == 0)
{
state._emat.remove_me(me, state._bg);
if (state._coupled_state != nullptr)
state._coupled_state->remove_edge(me);
else
boost::remove_edge(me, state._bg);
me = state._emat.get_null_edge();
}
});
};
if (state._rec_types.empty())
{
eops([](auto&&... args) { entries_op(args...);},
[](auto&){}, [](auto&){},
[](auto delta) { return delta == 0; });
if (state._coupled_state != nullptr)
{
m_entries._p_entries.clear();
std::vector<double> dummy;
entries_op(m_entries, state._emat,
[&](auto r, auto s, auto& me, auto delta)
{
if (delta == 0)
return;
m_entries._p_entries.emplace_back(r, s, me, delta,
dummy);
});
}
if (state._coupled_state != nullptr && !m_entries._p_entries.empty())
{
state._coupled_state->propagate_delta(m_entries.get_move().first,
m_entries.get_move().second,
m_entries._p_entries);
}
}
else
{
recs_apply_delta<Add, Remove>(state, m_entries, eops);
}
}
} // namespace graph_tool
#endif // GRAPH_BLOCKMODEL_ENTRIES_HH
src/graph/inference/blockmodel/graph_blockmodel_gibbs.hh
View file @ 4f0c7fe6
......@@ -95,11 +95,11 @@ struct Gibbs
{
if (nr == null_group)
{
if (!_allow_new_group)
if (!_allow_new_group || _state._allow_empty)
return numeric_limits<double>::infinity();
if (_state._empty_blocks.empty())
_state.add_block();
nr = _state._empty_blocks.front();
nr = _state._empty_blocks.back();
}
size_t r = _state._b[v];
if (!_state.allow_move(r, nr))
......@@ -107,13 +107,12 @@ struct Gibbs
return _state