Commit 1ade7acd authored by Tiago Peixoto's avatar Tiago Peixoto

Revamp multiflip mcmc and parametrize single node moves

parent 04dabc1c
Pipeline #287 failed with stage
in 238 minutes and 57 seconds
This diff is collapsed.
......@@ -117,8 +117,11 @@ for directed in [True, False]:
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=1, c=abs(c), niter=40)
if c < 0:
mcmc_args["w"] = g.vertex_index.copy()
else:
mcmc_args=dict(beta=1, niter=40)
if i == 0:
mcmc_equilibrate(state,
mcmc_args=mcmc_args,
......
......@@ -58,14 +58,11 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
double S = 0;
size_t nmoves = 0;
size_t nattempts = 0;
for (size_t iter = 0; iter < state._niter; ++iter)
{
if (!state._parallel)
{
std::shuffle(vlist.begin(), vlist.end(), rng_);
}
else
if (state._parallel)
{
parallel_loop(vlist,
[&](size_t, auto v)
......@@ -75,9 +72,14 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
numeric_limits<double>::max());
});
}
else
{
if (!state._deterministic)
std::shuffle(vlist.begin(), vlist.end(), rng_);
}
#pragma omp parallel firstprivate(state, probs, deltas, idx) \
if (state._parallel)
reduction(+: S, nmoves, nattempts) if (state._parallel)
parallel_loop_no_spawn
(vlist,
[&](size_t, auto v)
......@@ -92,6 +94,8 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
auto& moves = state.get_moves(v);
nattempts += moves.size();
probs.resize(moves.size());
deltas.resize(moves.size());
idx.resize(moves.size());
......@@ -167,7 +171,7 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
}
}
}
return make_pair(S, nmoves);
return std::make_tuple(S, nmoves, nattempts);
}
} // graph_tool namespace
......
......@@ -1152,21 +1152,21 @@ public:
// Sample node placement
template <class RNG>
size_t sample_block(size_t v, double c, RNG& rng)
size_t sample_block(size_t v, double c, double d, RNG& rng)
{
std::uniform_real_distribution<> rdist;
// attempt random block
size_t s;
if (_empty_blocks.empty())
if (!_empty_blocks.empty() && rdist(rng) < d)
{
s = uniform_sample(_candidate_blocks.begin() + 1,
_candidate_blocks.end(),
rng);
return uniform_sample(_empty_blocks, rng);
}
else
{
s = uniform_sample(_candidate_blocks, rng);
if (s == null_group)
s = uniform_sample(_empty_blocks, rng);
s = uniform_sample(_candidate_blocks.begin() + 1,
_candidate_blocks.end(),
rng);
}
if (!std::isinf(c) && !_neighbour_sampler.empty(v))
......@@ -1176,16 +1176,14 @@ public:
double p_rand = 0;
if (c > 0)
{
size_t B = (_empty_blocks.empty()) ?
_candidate_blocks.size() - 1 : _candidate_blocks.size();
size_t B = _candidate_blocks.size() - 1;
if (is_directed::apply<g_t>::type::value)
p_rand = c * B / double(_mrp[t] + _mrm[t] + c * B);
else
p_rand = c * B / double(_mrp[t] + c * B);
}
typedef std::uniform_real_distribution<> rdist_t;
if (c == 0 || rdist_t()(rng) >= p_rand)
if (c == 0 || rdist(rng) >= p_rand)
{
if (_egroups.empty())
_egroups.init(_b, _eweight, _g, _bg);
......@@ -1201,9 +1199,9 @@ public:
return s;
}
size_t sample_block(size_t v, double c, rng_t& rng)
size_t sample_block(size_t v, double c, double d, rng_t& rng)
{
return sample_block<rng_t>(v, c, rng);
return sample_block<rng_t>(v, c, d, rng);
}
size_t random_neighbour(size_t v, rng_t& rng)
......@@ -1215,12 +1213,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, bool reverse,
MEntries& m_entries)
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse, MEntries& m_entries)
{
typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
size_t B = (_empty_blocks.empty()) ?
_candidate_blocks.size() - 1 : _candidate_blocks.size();
size_t B = _candidate_blocks.size() - 1;
if (reverse)
{
if (_wr[s] == _vweight[v])
return d;
if (_wr[r] == 0)
B++;
// if (_wr[s] == _vweight[v])
// B--;
}
else
{
if (_wr[s] == 0)
return d;
}
if (B == num_vertices(_bg))
d = 0;
if (std::isinf(c))
return (1. - d) / B;
double p = 0;
size_t w = 0;
......@@ -1232,7 +1250,7 @@ public:
auto sum_prob = [&](auto& e, auto u)
{
vertex_t t = _b[u];
size_t t = _b[u];
if (u == v)
t = r;
size_t ew = _eweight[e];
......@@ -1306,15 +1324,16 @@ 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, bool reverse)
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse)
{
get_move_entries(v, _b[v], (reverse) ? r : s, _m_entries);
return get_move_prob(v, r, s, c, reverse, _m_entries);
return get_move_prob(v, r, s, c, d, reverse, _m_entries);
}
bool is_last(size_t v)
......
......@@ -49,7 +49,8 @@ python::object do_gibbs_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
});
};
block_state::dispatch(oblock_state, dispatch);
......
......@@ -42,6 +42,7 @@ using namespace std;
((allow_vacate,, bool, 0)) \
((parallel,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
((niter,, size_t, 0))
......
......@@ -52,12 +52,12 @@ void export_sbm_state()
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, rng_t&) =
size_t (state_t::*sample_block)(size_t, double, double, rng_t&) =
&state_t::sample_block;
size_t (state_t::*random_neighbour)(size_t, rng_t&) =
&state_t::random_neighbour;
double (state_t::*get_move_prob)(size_t, size_t, size_t, double,
bool) =
double, bool) =
&state_t::get_move_prob;
void (state_t::*merge_vertices)(size_t, size_t) =
&state_t::merge_vertices;
......
......@@ -421,7 +421,7 @@ struct Layers
}
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c,
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse, MEntries& m_entries)
{
// m_entries may include entries from different levels
......@@ -430,13 +430,13 @@ struct Layers
m_entries.clear();
BaseState::get_move_entries(v, r, s, m_entries);
}
return BaseState::get_move_prob(v, r, s, c, reverse, m_entries);
return BaseState::get_move_prob(v, r, s, c, d, reverse, m_entries);
}
double get_move_prob(size_t v, size_t r, size_t s, double c,
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse)
{
return BaseState::get_move_prob(v, r, s, c, reverse);
return BaseState::get_move_prob(v, r, s, c, d, reverse);
}
void merge_vertices(size_t u, size_t v)
......
......@@ -62,7 +62,8 @@ python::object gibbs_layered_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
});
},
false);
......
......@@ -53,10 +53,10 @@ void export_lsbm()
double (state_t::*virtual_move)(size_t, size_t, size_t, entropy_args_t) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, rng_t&)
size_t (state_t::*sample_block)(size_t, double, double, rng_t&)
= &state_t::sample_block;
double (state_t::*get_move_prob)(size_t, size_t, size_t, double,
bool)
double, bool)
= &state_t::get_move_prob;
void (state_t::*merge_vertices)(size_t, size_t)
= &state_t::merge_vertices;
......
......@@ -70,17 +70,20 @@ void export_layered_overlap_blockmodel_state()
layered_block_state<block_state_t>::dispatch
([&](auto* s)
{
typedef typename std::remove_reference<decltype(*s)>::type state_t;
typedef typename std::remove_reference<decltype(*s)>::type
state_t;
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, rng_t&)
size_t (state_t::*sample_block)(size_t, double, double,
rng_t&)
= &state_t::sample_block;
double (state_t::*get_move_prob)(size_t, size_t, size_t, double,
bool)
double (state_t::*get_move_prob)(size_t, size_t, size_t,
double, double, bool)
= &state_t::get_move_prob;
void (state_t::*move_vertices)(python::object, python::object) =
void (state_t::*move_vertices)(python::object,
python::object) =
&state_t::move_vertices;
class_<state_t> c(name_demangle(typeid(state_t).name()).c_str(),
......
......@@ -62,7 +62,8 @@ python::object gibbs_layered_overlap_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
});
},
false);
......
......@@ -39,10 +39,12 @@ using namespace std;
((vlist,&, std::vector<size_t>&, 0)) \
((beta,, double, 0)) \
((c,, double, 0)) \
((d,, double, 0)) \
((entropy_args,, entropy_args_t, 0)) \
((allow_vacate,, bool, 0)) \
((parallel,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
((niter,, size_t, 0))
......@@ -83,6 +85,11 @@ struct MCMC
return _state._b[v];
}
bool skip_node(size_t v)
{
return _state.node_weight(v) == 0;
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
......@@ -94,12 +101,12 @@ struct MCMC
auto r = _state._b[v];
if (!_allow_vacate && _state.is_last(v))
return r;
return null_group;
size_t s = _state.sample_block(v, _c, rng);
size_t s = _state.sample_block(v, _c, _d, rng);
if (!_state.allow_move(r, s))
return r;
return null_group;
return s;
}
......@@ -111,11 +118,11 @@ struct MCMC
double dS = _state.virtual_move(v, r, nr, _entropy_args,
_m_entries);
double a = 0;
if (!std::isinf(_c) && !std::isinf(_beta))
if (!std::isinf(_beta))
{
double pf = _state.get_move_prob(v, r, nr, _c, false,
double pf = _state.get_move_prob(v, r, nr, _c, _d, false,
_m_entries);
double pb = _state.get_move_prob(v, nr, r, _c, true,
double pb = _state.get_move_prob(v, nr, r, _c, _d, true,
_m_entries);
a = log(pb) - log(pf);
}
......
......@@ -46,10 +46,11 @@ using namespace std;
((E,, size_t, 0)) \
((vlist,&, std::vector<size_t>&, 0)) \
((c,, double, 0)) \
((d,, double, 0)) \
((entropy_args,, entropy_args_t, 0)) \
((allow_vacate,, bool, 0)) \
((verbose,, bool, 0)) \
((target_bin,, int, 0)) \
((target_bin,, int, 0)) \
((niter,, size_t, 0))
......@@ -104,7 +105,7 @@ struct Multicanonical
{
auto r = _state._b[v];
size_t s = _state.sample_block(v, _c, rng);
size_t s = _state.sample_block(v, _c, _d, rng);
if (!_state.allow_move(r, s))
return r;
......@@ -119,14 +120,10 @@ struct Multicanonical
{
double dS = _state.virtual_move(v, _state._b[v], nr, _entropy_args);
double a = 0;
if (!std::isinf(_c))
{
size_t r = _state._b[v];
double pf = _state.get_move_prob(v, r, nr, _c, false);
double pb = _state.get_move_prob(v, nr, r, _c, true);
a = log(pb) - log(pf);
}
size_t r = _state._b[v];
double pf = _state.get_move_prob(v, r, nr, _c, _d, false);
double pb = _state.get_move_prob(v, nr, r, _c, _d, true);
double a = log(pb) - log(pf);
return make_pair(dS, a);
}
......
......@@ -56,8 +56,62 @@ python::object do_multiflip_mcmc_sweep(python::object omcmc_state,
return ret;
}
void get_worder(GraphInterface& gi, boost::any abs, boost::any aorder)
{
typedef vprop_map_t<int64_t>::type vmap_t;
vmap_t& order = any_cast<vmap_t&>(aorder);
run_action<>()(gi,
[&](auto& g, auto& bs)
{
std::vector<size_t> vs;
vs.reserve(num_vertices(g));
for (auto v : vertices_range(g))
vs.push_back(v);
std::sort(vs.begin(), vs.end(),
[&](auto u, auto v)
{ return bs[u] < bs[v]; });
for (size_t i = 0; i < vs.size(); ++i)
order[vs[i]] = i;
},
vertex_scalar_vector_properties())(abs);
}
void shuffle_labels(GraphInterface& gi, boost::any ab, rng_t& rng)
{
run_action<>()(gi,
[&](auto& g, auto& b)
{
typedef typename property_traits
<typename std::remove_reference<decltype(b)>::type>
::value_type val_t;
std::unordered_set<val_t> vals;
for (auto v : vertices_range(g))
vals.insert(b[v]);
std::vector<val_t> nvals(vals.begin(), vals.end());
std::shuffle(nvals.begin(), nvals.end(), rng);
std::unordered_map<val_t, val_t> val_map;
auto iter = nvals.begin();
for (auto r : vals)
val_map[r] = *iter++;
for (auto v : vertices_range(g))
b[v] = val_map[b[v]];
},
writable_vertex_scalar_properties())(ab);
}
void export_blockmodel_multiflip_mcmc()
{
using namespace boost::python;
def("multiflip_mcmc_sweep", &do_multiflip_mcmc_sweep);
def("get_worder", &get_worder);
def("shuffle_labels", &shuffle_labels);
}
......@@ -158,10 +158,10 @@ void export_overlap_blockmodel_state()
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, rng_t&)
size_t (state_t::*sample_block)(size_t, double, double, rng_t&)
= &state_t::sample_block;
double (state_t::*get_move_prob)(size_t, size_t, size_t, double,
bool)
double, bool)
= &state_t::get_move_prob;
void (state_t::*set_partition)(boost::any&)
= &state_t::set_partition;
......
......@@ -557,7 +557,7 @@ public:
// Sample node placement
template <class RNG>
size_t sample_block(size_t v, double c, RNG& rng)
size_t sample_block(size_t v, double c, double, RNG& rng)
{
// attempt random block
size_t s = uniform_sample(_candidate_blocks, rng);
......@@ -596,9 +596,9 @@ public:
return s;
}
size_t sample_block(size_t v, double c, rng_t& rng)
size_t sample_block(size_t v, double c, double d, rng_t& rng)
{
return sample_block<rng_t>(v, c, rng);
return sample_block<rng_t>(v, c, d, rng);
}
template <class RNG>
......@@ -622,8 +622,8 @@ public:
// Computes the move proposal probability
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c, bool reverse,
MEntries& m_entries)
double get_move_prob(size_t v, size_t r, size_t s, double c, double,
bool reverse, MEntries& m_entries)
{
typedef typename graph_traits<g_t>::vertex_descriptor vertex_t;
size_t B = num_vertices(_bg);
......@@ -709,9 +709,10 @@ public:
return 1. / B;
}
double get_move_prob(size_t v, size_t r, size_t s, double c, bool reverse)
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse)
{
return get_move_prob(v, r, s, c, reverse, _m_entries);
return get_move_prob(v, r, s, c, d, reverse, _m_entries);
}
bool is_last(size_t v)
......
......@@ -53,7 +53,8 @@ python::object gibbs_overlap_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -39,9 +39,11 @@ using namespace std;
((vlist,, std::vector<size_t>, 0)) \
((beta,, double, 0)) \
((c,, double, 0)) \
((d,, double, 0)) \
((entropy_args,, entropy_args_t, 0)) \
((allow_vacate,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
((niter,, size_t, 0))
......@@ -110,6 +112,11 @@ struct MCMC
return _state._b[v];
}
size_t skip_node(size_t v)
{
return _state.node_weight(v) == 0;
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
......@@ -121,10 +128,10 @@ struct MCMC
auto v = _bundles[i][0];
auto r = _state._b[v];
size_t s = _state.sample_block(v, _c, rng);
size_t s = _state.sample_block(v, _c, _d, rng);
if (_state._bclabel[s] != _state._bclabel[r])
return r;
return null_group;
return s;
}
......
......@@ -121,7 +121,7 @@ struct Merge
else
{
auto v = uniform_sample(bundle, rng);
s = _state.sample_block(v, 0, rng);
s = _state.sample_block(v, 0, 0, rng);
}
if (s == r || _state._bclabel[r] != _state._bclabel[s])
......
......@@ -69,7 +69,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
for (size_t iter = 0; iter < state._niter; ++iter)
{
if (state._sequential)
if (state._sequential && !state._deterministic)
std::shuffle(vlist.begin(), vlist.end(), rng);
for (size_t vi = 0; vi < vlist.size(); ++vi)
......@@ -80,13 +80,13 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
else
v = uniform_sample(vlist, rng);
if (state.node_weight(v) == 0)
if (state.skip_node(v))
continue;