Commit c5e1c6a2 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

mcmc_equilibrate(): add "block_moves" argument

parent 1a33662e
Pipeline #134 failed with stage
......@@ -95,7 +95,7 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
v = vertex(vlist[v_rand(rng)], g);
}
if (state.skip(v))
if (state.node_weight(v) == 0)
continue;
vector<size_t>& moves = state.get_moves(v);
......@@ -146,7 +146,7 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
if (!state._parallel)
{
state.perform_move(v, s);
nmoves++;
nmoves += state.node_weight(v);
S += deltas[j];
}
else
......
......@@ -231,11 +231,15 @@ void export_blockmodel_state()
void (state_t::*merge_vertices)(size_t, size_t)
= &state_t::merge_vertices;
void (state_t::*set_partition)(boost::any&)
= &state_t::set_partition;
class_<state_t> c(name_demangle(typeid(state_t).name()).c_str(),
no_init);
c.def("remove_vertex", &state_t::remove_vertex)
.def("add_vertex", &state_t::add_vertex)
.def("move_vertex", &state_t::move_vertex)
.def("set_partition", set_partition)
.def("virtual_move", virtual_move)
.def("merge_vertices", merge_vertices)
.def("sample_block", sample_block)
......
......@@ -275,6 +275,19 @@ public:
add_vertex(v, nr);
}
template <class VMap>
void set_partition(VMap&& b)
{
for (auto v : vertices_range(_g))
move_vertex(v, b[v]);
}
void set_partition(boost::any& ab)
{
vmap_t& b = boost::any_cast<vmap_t&>(ab);
set_partition<typename vmap_t::unchecked_t>(b.get_unchecked());
}
size_t virtual_remove_size(size_t v)
{
return _wr[_b[v]] - _vweight[v];
......@@ -784,6 +797,11 @@ public:
return _wr[_b[v]] == _vweight[v];
}
size_t node_weight(size_t v)
{
return _vweight[v];
}
double get_deg_entropy(size_t v, const simple_degs_t&)
{
if (_ignore_degree[v] == 1)
......
......@@ -135,9 +135,10 @@ struct Gibbs
return _state._b[v];
}
bool skip(size_t v)
size_t node_weight(size_t v)
{
return _state._vweight[v] == 0;
return _state.node_weight(v);
}
double virtual_move_dS(size_t v, size_t nr)
......
......@@ -453,11 +453,15 @@ void export_layered_blockmodel_state()
void (state_t::*merge_vertices)(size_t, size_t)
= &state_t::merge_vertices;
void (state_t::*set_partition)(boost::any&)
= &state_t::set_partition;
class_<state_t> c(name_demangle(typeid(state_t).name()).c_str(),
no_init);
c.def("remove_vertex", &state_t::remove_vertex)
.def("add_vertex", &state_t::add_vertex)
.def("move_vertex", &state_t::move_vertex)
.def("set_partition", set_partition)
.def("virtual_move", virtual_move)
.def("merge_vertices", merge_vertices)
.def("sample_block", sample_block)
......
......@@ -217,6 +217,20 @@ struct Layers
}
}
template <class VMap>
void set_partition(VMap&& b)
{
for (auto v : vertices_range(_g))
LayeredBlockState::move_vertex(v, b[v]);
}
void set_partition(boost::any& ab)
{
typename BaseState::b_t::checked_t& b
= boost::any_cast<typename BaseState::b_t::checked_t&>(ab);
set_partition(b.get_unchecked());
}
template <class MEntries>
double virtual_move(size_t v, size_t s, bool dense, bool multigraph,
bool partition_dl, bool deg_dl, bool edges_dl,
......
......@@ -85,6 +85,11 @@ struct MCMC
return _state._b[v];
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
size_t move_proposal(size_t v, RNG& rng)
{
......
......@@ -95,6 +95,11 @@ struct Merge
return _state._b[v];
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
size_t move_proposal(size_t v, bool random, RNG& rng)
{
......
......@@ -87,6 +87,11 @@ struct Multicanonical
return _state._b[v];
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
size_t move_proposal(size_t v, RNG& rng)
{
......
......@@ -160,11 +160,15 @@ void export_overlap_blockmodel_state()
bool)
= &state_t::get_move_prob;
void (state_t::*set_partition)(boost::any&)
= &state_t::set_partition;
class_<state_t> c(name_demangle(typeid(state_t).name()).c_str(),
no_init);
c.def("remove_vertex", &state_t::remove_vertex)
.def("add_vertex", &state_t::add_vertex)
.def("move_vertex", &state_t::move_vertex)
.def("set_partition", set_partition)
.def("virtual_move", virtual_move)
.def("sample_block", sample_block)
.def("entropy", &state_t::entropy)
......
......@@ -217,6 +217,19 @@ public:
add_vertex(v, nr);
}
template <class VMap>
void set_partition(VMap&& b)
{
for (auto v : vertices_range(_g))
move_vertex(v, b[v]);
}
void set_partition(boost::any& ab)
{
vmap_t& b = boost::any_cast<vmap_t&>(ab);
set_partition<typename vmap_t::unchecked_t>(b.get_unchecked());
}
size_t virtual_remove_size(size_t v)
{
return _overlap_stats.virtual_remove_size(v, _b[v]);
......@@ -488,6 +501,11 @@ public:
return _overlap_stats.virtual_remove_size(v, r) == 0;
}
size_t node_weight(size_t)
{
return 1;
}
double sparse_entropy(bool multigraph, bool deg_entropy) const
{
double S = 0;
......
......@@ -113,6 +113,11 @@ struct MCMC
return _state._b[v];
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
size_t move_proposal(size_t i, RNG& rng)
{
......
......@@ -92,6 +92,9 @@ auto mcmc_sweep(MCMCState state, RNG& rng_)
v = vertex(vlist[v_rand(rng)], g);
}
if (state.node_weight(v) == 0)
continue;
auto r = state.node_state(v);
auto s = state.move_proposal(v, rng);
......@@ -125,7 +128,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng_)
if (!state._parallel)
{
state.perform_move(v, s);
nmoves++;
nmoves += state.node_weight(v);
S += dS.first;
}
else
......
......@@ -62,6 +62,9 @@ auto merge_sweep(MergeState state, RNG& rng_)
auto v = state._available[i];
if (state.node_weight(v) == 0)
continue;
gt_hash_set<size_t> past_moves;
auto find_candidates = [&](bool random)
......
......@@ -62,6 +62,9 @@ auto multicanonical_sweep(MulticanonicalState& state, RNG& rng)
{
auto v = vertex(uniform_sample(vlist, rng), g);
if (state.node_weight(v) == 0)
continue;
auto s = state.move_proposal(v, rng);
std::pair<double, double> dS = state.virtual_move_dS(v, s);
......
......@@ -346,6 +346,12 @@ class BlockState(object):
r"""Returns the property map which contains the block labels for each vertex."""
return self.b
def set_blocks(self, b):
r"""Sets the internal partition of the state."""
if b.value_type() != "int32_t":
b = b.copy("int32_t")
self._state.set_partition(_prop("v", self.g, b))
def get_bg(self):
r"""Returns the block graph."""
return self.bg
......@@ -732,7 +738,14 @@ class BlockState(object):
mcmc_state.block_list.extend(block_list)
mcmc_state.vlist = Vector_size_t()
if vertices is None:
mcmc_state.vlist.extend(self.g.vertex_index.copy().fa)
idx = self.g.vertex_index.copy().fa
if (hasattr(self, "vweight") and
not isinstance(self.vweight, libinference.unity_vprop_t)):
# ignore vertices with zero weight
vw = self.vweight.fa
mcmc_state.vlist.extend(idx[vw > 0])
else:
mcmc_state.vlist.extend(idx)
else:
mcmc_state.vlist.extend(vertices)
mcmc_state.E = self.E
......@@ -821,7 +834,14 @@ class BlockState(object):
gibbs_state.block_list.extend(block_list)
gibbs_state.vlist = Vector_size_t()
if vertices is None:
gibbs_state.vlist.extend(self.g.vertex_index.copy().fa)
idx = self.g.vertex_index.copy().fa
if (hasattr(self, "vweight") and
not isinstance(self.vweight, libinference.unity_vprop_t)):
# ignore vertices with zero weight
vw = self.vweight.fa
gibbs_state.vlist.extend(idx[vw > 0])
else:
gibbs_state.vlist.extend(idx)
else:
gibbs_state.vlist.extend(vertices)
gibbs_state.E = self.E
......@@ -915,7 +935,14 @@ class BlockState(object):
multi_state.block_list.extend(block_list)
multi_state.vlist = Vector_size_t()
if vertices is None:
multi_state.vlist.extend(self.g.vertex_index.copy().fa)
idx = self.g.vertex_index.copy().fa
if (hasattr(self, "vweight") and
not isinstance(self.vweight, libinference.unity_vprop_t)):
# ignore vertices with zero weight
vw = self.vweight.fa
multi_state.vlist.extend(idx[vw > 0])
else:
multi_state.vlist.extend(idx)
else:
multi_state.vlist.extend(vertices)
multi_state.E = self.E
......
......@@ -29,9 +29,9 @@ import numpy
from . util import *
def mcmc_equilibrate(state, wait=10, nbreaks=2, max_niter=numpy.inf,
force_niter=None, epsilon=0, gibbs=False, mcmc_args={},
entropy_args={}, history=False, callback=None,
verbose=False):
force_niter=None, epsilon=0, gibbs=False,
block_moves=False, mcmc_args={}, entropy_args={},
history=False, callback=None, verbose=False):
r"""Equilibrate a MCMC with a given starting state.
Parameters
......@@ -54,6 +54,9 @@ def mcmc_equilibrate(state, wait=10, nbreaks=2, max_niter=numpy.inf,
gibbs : ``bool`` (optional, default: ``False``)
If ``True``, each step will call ``state.gibbs_sweep`` instead of
``state.mcmc_sweep``.
block_moves : ``bool`` (optional, default: ``False``)
If ``True``, each iteration will be accompanied by a "block move", where
all vertices of the same group are moved simultaneously.
mcmc_args : ``dict`` (optional, default: ``{}``)
Arguments to be passed to ``state.mcmc_sweep`` (or ``state.gibbs_sweep``).
history : ``bool`` (optional, default: ``False``)
......@@ -107,6 +110,21 @@ def mcmc_equilibrate(state, wait=10, nbreaks=2, max_niter=numpy.inf,
else:
delta, nmoves = state.gibbs_sweep(**mcmc_args)
if block_moves:
bstate = state.get_block_state(vweight=True,
clabel=state.get_bclabel())
if not gibbs:
ret = bstate.mcmc_sweep(**mcmc_args)
else:
ret = bstate.gibbs_sweep(**mcmc_args)
b = state.b.copy()
pmap(b, bstate.b)
state.set_blocks(b)
delta += ret[0]
nmoves += ret[1]
S += delta
niter += 1
total_nmoves += 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