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

General improvements to multiflip and gibbs mcmc

parent 9b24872c
Pipeline #529 failed with stage
in 619 minutes and 5 seconds
......@@ -117,16 +117,13 @@ for nested in [False, True]:
state = minimize_blockmodel_dl(g, deg_corr=True)
state = state.copy(B=g.num_vertices())
if nested:
cs = list(reversed([-numpy.inf, numpy.inf, 1, 0.1, 0.01, -1]))
else:
cs = list(reversed([-numpy.inf, numpy.inf, 1, 0.1, 0.01, "gibbs", -1]))
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs", -numpy.inf, -1]))
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=.2, c=abs(c), niter=40)
mcmc_args=dict(beta=.3, c=abs(c), niter=40)
else:
mcmc_args=dict(beta=.2, niter=40)
mcmc_args=dict(beta=.3, niter=40)
if nested:
get_B = lambda s: [sl.get_nonempty_B() for sl in s.levels]
......@@ -146,8 +143,8 @@ for nested in [False, True]:
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
wait=1000,
nbreaks=5,
wait=6000,
nbreaks=10,
callback=get_B,
verbose=(1, "c = %s " % str(c)) if verbose else False,
history=True)
......
......@@ -68,53 +68,53 @@ public:
typedef checked_vector_property_map<T,IndexMap> self_t;
checked_vector_property_map(const IndexMap& idx = IndexMap())
: store(std::make_shared<std::vector<T>>()), index(idx) {}
: _store(std::make_shared<std::vector<T>>()), _index(idx) {}
checked_vector_property_map(unsigned initial_size,
const IndexMap& idx = IndexMap())
: store(std::make_shared<std::vector<T>>(initial_size)), index(idx) {}
: _store(std::make_shared<std::vector<T>>(initial_size)), _index(idx) {}
typename std::vector<T>::iterator storage_begin()
{
return store->begin();
return _store->begin();
}
typename std::vector<T>::iterator storage_end()
{
return store->end();
return _store->end();
}
typename std::vector<T>::const_iterator storage_begin() const
{
return store->begin();
return _store->begin();
}
typename std::vector<T>::const_iterator storage_end() const
{
return store->end();
return _store->end();
}
void reserve(size_t size) const
{
if (size > store->size())
store->resize(size);
if (size > _store->size())
_store->resize(size);
}
void resize(size_t size) const
{
store->resize(size);
_store->resize(size);
}
void shrink_to_fit() const
{
store->shrink_to_fit();
_store->shrink_to_fit();
}
std::vector<T>& get_storage() const { return (*store); }
std::vector<T>& get_storage() const { return *_store; }
void swap(checked_vector_property_map& other)
{
store->swap(*other.store);
_store->swap(*other._store);
}
unchecked_t get_unchecked(size_t size = 0) const
......@@ -126,32 +126,23 @@ public:
// deep copy
checked_vector_property_map copy() const
{
checked_vector_property_map pmap(index);
*(pmap.store) = *store;
checked_vector_property_map pmap(_index);
*(pmap._store) = *_store;
return pmap;
}
public:
// Copy ctor absent, default semantics is OK.
// Assignment operator absent, default semantics is OK.
// CONSIDER: not sure that assignment to 'index' is correct.
reference operator[](const key_type& v) const {
auto i = get(index, v);
if (static_cast<size_t>(i) >= store->size()) {
store->resize(i + 1, T());
auto i = get(_index, v);
auto& store = *_store;
if (static_cast<size_t>(i) >= store.size()) {
store.resize(i + 1);
}
return (*store)[i];
return store[i];
}
protected:
// Conceptually, we have a vector of infinite size. For practical
// purposes, we start with an empty vector and grow it as needed.
// Note that we cannot store pointer to vector here -- we cannot
// store pointer to data, because if copy of property map resizes
// the vector, the pointer to data will be invalidated.
// I wonder if class 'pmap_ref' is simply needed.
std::shared_ptr< std::vector<T> > store;
IndexMap index;
std::shared_ptr<std::vector<T>> _store;
IndexMap _index;
};
template<typename T, typename IndexMap = identity_property_map>
......@@ -174,14 +165,14 @@ public:
size_t size = 0)
: _checked(checked)
{
if (size > 0 && _checked.store->size() < size)
_checked.store->resize(size);
if (size > 0 && _checked._store->size() < size)
_checked._store->resize(size);
}
unchecked_vector_property_map(const IndexMap& index_map,
size_t size = 0)
: _checked(size, index_map)
{
*this = unchecked_vector_property_map(checked_t(index_map), size);
}
void reserve(size_t size) const { _checked.reserve(size); }
......@@ -192,7 +183,7 @@ public:
__attribute__((always_inline)) __attribute__((flatten))
reference operator[](const key_type& v) const
{
return (*_checked.store)[get(_checked.index, v)];
return (*_checked._store)[get(_checked._index, v)];
}
std::vector<T>& get_storage() const { return _checked.get_storage(); }
......
......@@ -192,6 +192,10 @@ public:
}
_dBdx.resize(_rec_types.size());
_LdBdx.resize(_rec_types.size());
_N = 0;
for (auto v : vertices_range(_g))
_N += _vweight[v];
}
BlockState(const BlockState& other)
......@@ -209,6 +213,7 @@ public:
_B_E(other._B_E),
_B_E_D(other._B_E_D),
_rt(other._rt),
_N(other._N),
_vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
_eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
_degs(uncheck(__adegs, typename std::add_pointer<degs_t>::type())),
......@@ -275,15 +280,16 @@ public:
BlockState::remove_partition_node(v, r);
}
bool allow_move(size_t v, size_t r, size_t nr)
bool allow_move(size_t r, size_t nr)
{
if (_coupled_state != nullptr && is_last(v))
if (_coupled_state != nullptr)
{
auto& bh = _coupled_state->get_b();
if (bh[r] != bh[nr])
auto& hb = _coupled_state->get_b();
auto rr = hb[r];
auto ss = hb[nr];
if (rr != ss && !_coupled_state->allow_move(rr, ss))
return false;
}
return _bclabel[r] == _bclabel[nr];
}
......@@ -293,7 +299,7 @@ public:
if (r == nr)
return;
if (!allow_move(v, r, nr))
if (!allow_move(r, nr))
throw ValueException("cannot move vertex across clabel barriers");
get_move_entries(v, r, nr, _m_entries, std::forward<EFilt>(efilt));
......@@ -861,7 +867,9 @@ public:
template <class VMap>
void set_vertex_weight(size_t v, int w, VMap&& vweight)
{
_N -= vweight[v];
vweight[v] = w;
_N += w;
}
void init_vertex_weight(size_t v)
......@@ -1384,7 +1392,7 @@ public:
{
assert(size_t(_b[v]) == r || r == null_group);
if (r != null_group && nr != null_group && !allow_move(v, r, nr))
if (r != null_group && nr != null_group && !allow_move(r, nr))
return std::numeric_limits<double>::infinity();
get_move_entries(v, r, nr, m_entries, [](auto) { return false; });
......@@ -1722,17 +1730,12 @@ public:
// Move proposals
// =========================================================================
// Sample node placement
size_t sample_block(size_t v, double c, double d, rng_t& rng)
size_t get_empty_block(size_t v)
{
// attempt new block
size_t s;
std::bernoulli_distribution new_r(d);
if (d > 0 && new_r(rng) && (_candidate_blocks.size() - 1 < num_vertices(_g)))
if (_empty_blocks.empty())
{
if (_empty_blocks.empty())
add_block();
s = uniform_sample(_empty_blocks, rng);
add_block();
auto s = _empty_blocks.back();
auto r = _b[v];
_bclabel[s] = _bclabel[r];
if (_coupled_state != nullptr)
......@@ -1740,9 +1743,51 @@ public:
auto& hb = _coupled_state->get_b();
hb[s] = hb[r];
}
}
return _empty_blocks.back();
}
void sample_branch(size_t v, size_t u, rng_t& rng)
{
size_t s;
std::bernoulli_distribution new_r(1./_candidate_blocks.size());
if (_candidate_blocks.size() <= num_vertices(_g) && new_r(rng))
{
get_empty_block(v);
s = uniform_sample(_empty_blocks, rng);
auto r = _b[u];
if (_coupled_state != nullptr)
_coupled_state->sample_branch(s, r, rng);
_bclabel[s] = _bclabel[r];
}
else
{
s = uniform_sample(_candidate_blocks.begin() + 1,
_candidate_blocks.end(), rng);
}
_b[v] = s;
}
// Sample node placement
size_t sample_block(size_t v, double c, double d, rng_t& rng)
{
size_t B = _candidate_blocks.size() - 1;
// attempt new block
std::bernoulli_distribution new_r(d);
if (d > 0 && B < _N && new_r(rng))
{
get_empty_block(v);
auto s = uniform_sample(_empty_blocks, rng);
auto r = _b[v];
if (_coupled_state != nullptr)
_coupled_state->sample_branch(s, r, rng);
_bclabel[s] = _bclabel[r];
return s;
}
size_t s;
if (!std::isinf(c) && !_neighbor_sampler.empty(v))
{
auto u = _neighbor_sampler.sample(v, rng);
......@@ -1750,7 +1795,6 @@ public:
double p_rand = 0;
if (c > 0)
{
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
......@@ -1791,6 +1835,7 @@ public:
return _neighbor_sampler.sample(v, rng);
}
// Computes the move proposal probability
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
......@@ -1804,8 +1849,6 @@ public:
return d;
if (_wr[r] == 0)
B++;
// if (_wr[s] == _vweight[v])
// B--;
}
else
{
......@@ -1813,7 +1856,7 @@ public:
return d;
}
if (B == num_vertices(_g))
if (B == _N)
d = 0;
if (std::isinf(c))
......@@ -1824,11 +1867,11 @@ public:
size_t kout = 0, kin = 0;
degs_op(v, _vweight, _eweight, _degs, _g,
[&] ([[maybe_unused]] size_t din, size_t dout, auto c)
[&] ([[maybe_unused]] size_t din, size_t dout, auto count)
{
kout += dout * c;
kout += dout * count;
if constexpr (is_directed_::apply<g_t>::type::value)
kin += din * c;
kin += din * count;
});
if constexpr (!is_directed_::apply<g_t>::type::value)
kin = kout;
......@@ -1939,7 +1982,7 @@ public:
bool is_last(size_t v)
{
return _wr[_b[v]] == _vweight[v];
return (_vweight[v] > 0) && (_wr[_b[v]] == _vweight[v]);
}
size_t node_weight(size_t v)
......@@ -2570,6 +2613,7 @@ public:
size_t _B_E = 0;
size_t _B_E_D = 0;
int _rt = weight_type::NONE;
size_t _N;
typedef typename std::conditional<is_weighted_t::value,
vmap_t::unchecked_t, vcmap_t>::type vweight_t;
......
......@@ -41,7 +41,6 @@ using namespace std;
((entropy_args,, entropy_args_t, 0)) \
((allow_vacate,, bool, 0)) \
((allow_new_group,, bool, 0)) \
((parallel,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
......@@ -67,7 +66,6 @@ struct Gibbs
sizeof...(Ts)>* = nullptr>
GibbsBlockState(ATs&&... as)
: GibbsBlockStateBase<Ts...>(as...),
_g(_state._g),
_m_entries(num_vertices(_state._bg))
{
_state.init_mcmc(numeric_limits<double>::infinity(),
......@@ -76,10 +74,14 @@ struct Gibbs
_entropy_args.edges_dl));
}
typename state_t::g_t& _g;
typename state_t::m_entries_t _m_entries;
auto& get_moves(size_t) { return _state._candidate_blocks; }
std::vector<size_t> _candidate_blocks;
auto& get_moves(size_t)
{
return _state._candidate_blocks;
}
size_t node_state(size_t v)
{
......@@ -91,32 +93,29 @@ struct Gibbs
return _state.node_weight(v);
}
double virtual_move_dS(size_t v, size_t nr)
size_t _nr;
double virtual_move_dS(size_t v, size_t nr, rng_t& rng)
{
size_t r = _state._b[v];
if (nr == null_group)
{
if (!_allow_new_group)
if (!_allow_new_group ||
_state._candidate_blocks.size() - 1 == num_vertices(_state._g))
return numeric_limits<double>::infinity();
if (_state._empty_blocks.empty())
_state.add_block();
nr = _state._empty_blocks.back();
_state.get_empty_block(v);
_nr = nr = uniform_sample(_state._empty_blocks, rng);
if (_state._coupled_state != nullptr)
_state._coupled_state->sample_branch(nr, r, rng);
_state._bclabel[nr] = _state._bclabel[r];
}
size_t r = _state._b[v];
if (!_state.allow_move(v, r, nr))
return numeric_limits<double>::infinity();
return _state.virtual_move(v, r, nr, _entropy_args, _m_entries);
}
void perform_move(size_t v, size_t nr)
{
size_t r = _state._b[v];
if (nr == null_group)
nr = _state._empty_blocks.back();
if (r == nr)
return;
nr = _nr;
_state.move_vertex(v, nr);
}
};
......
......@@ -48,16 +48,8 @@ python::object do_mcmc_sweep(python::object omcmc_state,
(omcmc_state,
[&](auto& s)
{
if (s._parallel)
{
auto ret_ = mcmc_sweep_parallel(s, rng);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
else
{
auto ret_ = mcmc_sweep(s, rng);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
auto ret_ = mcmc_sweep(s, rng);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
......
......@@ -42,7 +42,6 @@ using namespace std;
((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)) \
......@@ -79,7 +78,7 @@ struct MCMC
typename state_t::g_t& _g;
typename state_t::m_entries_t _m_entries;
size_t _null_move = null_group;
constexpr static size_t _null_move = null_group;
size_t node_state(size_t v)
{
......@@ -91,25 +90,12 @@ struct MCMC
return _state.node_weight(v) == 0;
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
size_t move_proposal(size_t v, RNG& rng)
{
auto r = _state._b[v];
if (!_allow_vacate && _state.is_last(v))
return null_group;
size_t s = _state.sample_block(v, _c, _d, rng);
if (!_state.allow_move(v, r, s))
return null_group;
return s;
return _null_move;
return _state.sample_block(v, _c, _d, rng);
}
std::tuple<double, double>
......
......@@ -113,7 +113,7 @@ struct Merge
s = uniform_sample(_available, rng);
}
if (s == r || !_state.allow_move(v, r, s))
if (s == r || !_state.allow_move(r, s))
return _null_move;
return s;
......
......@@ -91,13 +91,8 @@ struct Multicanonical
return _state.node_state(v);
}
size_t node_weight(size_t v)
{
return _state.node_weight(v);
}
template <class RNG>
move_t move_proposal(size_t v, RNG& rng)
auto move_proposal(size_t v, RNG& rng)
{
return _state.move_proposal(v, rng);
}
......@@ -141,6 +136,11 @@ struct Multicanonical
return _state.get_vlist();
}
size_t get_N()
{
return _state.get_N();
}
double get_beta()
{
return 1;
......
......@@ -92,14 +92,14 @@ struct MCMC
if (_state._vweight[v] == 0)
continue;
add_element(_groups[_state._b[v]], _vpos, v);
_N++;
_N += _state._vweight[v];
}
for (auto r : vertices_range(_state._bg))
{
if (_state._wr[r] == 0)
continue;
add_element(_vlist, _rpos, r);
add_element(_rlist, _rpos, r);
}
}
......@@ -114,13 +114,11 @@ struct MCMC
typename vprop_map_t<int>::type::unchecked_t _bprev;
typename vprop_map_t<int>::type::unchecked_t _bnext;
std::vector<size_t> _vlist;
std::vector<size_t> _rlist;
std::vector<size_t> _vs;
std::vector<size_t> _vs_temp;
constexpr static move_t _null_move = move_t::null;
size_t _nproposals = 0;
size_t _N = 0;
double _dS;
......@@ -135,20 +133,22 @@ struct MCMC
return r;
}
bool skip_node(size_t r)
constexpr bool skip_node(size_t)
{
return _groups[r].empty();
return false;
}