If you want to fork this repository to propose merge requests, please send an email to tiago@skewed.de, and your project limit will be raised.

Commit 5239342c authored by Tiago Peixoto's avatar Tiago Peixoto

inference: Add multiflip MCMC sweep

parent 93d7eec3
Pipeline #222 failed with stage
......@@ -112,23 +112,25 @@ for directed in [True, False]:
state = minimize_blockmodel_dl(g, deg_corr=True)
state = state.copy(B=g.num_vertices())
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs"]))
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs", -1]))
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=1, c=c, niter=40)
mcmc_args=dict(beta=1, c=abs(c), niter=40)
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,
......
......@@ -27,16 +27,19 @@ libgraph_tool_inference_la_SOURCES = \
graph_blockmodel_layers_mcmc.cc \
graph_blockmodel_layers_merge.cc \
graph_blockmodel_layers_multicanonical.cc \
graph_blockmodel_layers_multiflip_mcmc.cc \
graph_blockmodel_layers_overlap.cc \
graph_blockmodel_layers_overlap_exhaustive.cc \
graph_blockmodel_layers_overlap_mcmc.cc \
graph_blockmodel_layers_overlap_mcmc_bundled.cc \
graph_blockmodel_layers_overlap_multiflip_mcmc.cc \
graph_blockmodel_layers_overlap_gibbs.cc \
graph_blockmodel_layers_overlap_multicanonical.cc \
graph_blockmodel_layers_overlap_vacate.cc \
graph_blockmodel_marginals.cc \
graph_blockmodel_mcmc.cc \
graph_blockmodel_multicanonical.cc \
graph_blockmodel_multiflip_mcmc.cc \
graph_blockmodel_merge.cc \
graph_blockmodel_overlap.cc \
graph_blockmodel_overlap_exhaustive.cc \
......@@ -44,6 +47,7 @@ libgraph_tool_inference_la_SOURCES = \
graph_blockmodel_overlap_mcmc.cc \
graph_blockmodel_overlap_mcmc_bundled.cc \
graph_blockmodel_overlap_multicanonical.cc \
graph_blockmodel_overlap_multiflip_mcmc.cc \
graph_blockmodel_overlap_vacate.cc \
graph_inference.cc \
int_part.cc \
......@@ -63,6 +67,7 @@ libgraph_tool_inference_la_include_HEADERS = \
graph_blockmodel_mcmc.hh \
graph_blockmodel_merge.hh \
graph_blockmodel_multicanonical.hh \
graph_blockmodel_multiflip_mcmc.hh \
graph_blockmodel_overlap.hh \
graph_blockmodel_overlap_mcmc_bundled.hh \
graph_blockmodel_overlap_util.hh \
......
......@@ -114,6 +114,7 @@ public:
_vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
_eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
_emat(_bg, rng),
_egroups_enabled(true),
_neighbour_sampler(_g, _eweight),
_m_entries(num_vertices(_bg)),
_coupled_state(nullptr)
......@@ -139,6 +140,7 @@ public:
_vweight(uncheck(__avweight, typename std::add_pointer<vweight_t>::type())),
_eweight(uncheck(__aeweight, typename std::add_pointer<eweight_t>::type())),
_emat(other._emat),
_egroups_enabled(other._egroups_enabled),
_neighbour_sampler(other._neighbour_sampler),
_m_entries(num_vertices(_bg)),
_coupled_state(nullptr)
......@@ -147,7 +149,6 @@ public:
enable_partition_stats();
}
// =========================================================================
// State modification
// =========================================================================
......@@ -264,7 +265,7 @@ public:
{
_wr[r] -= _vweight[v];
if (!_egroups.empty())
if (!_egroups.empty() && _egroups_enabled)
_egroups.remove_vertex(v, _b, _g);
if (is_partition_stats_enabled())
......@@ -283,7 +284,7 @@ public:
{
_wr[r] += _vweight[v];
if (!_egroups.empty())
if (!_egroups.empty() && _egroups_enabled)
_egroups.add_vertex(v, _b, _eweight, _g);
if (is_partition_stats_enabled())
......@@ -1677,6 +1678,7 @@ public:
emat_t _emat;
EGroups<g_t, is_weighted_t> _egroups;
bool _egroups_enabled;
typedef NeighbourSampler<g_t, is_weighted_t, boost::mpl::false_>
neighbour_sampler_t;
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_blockmodel_util.hh"
#include "graph_blockmodel.hh"
#define BASE_STATE_params BLOCK_STATE_params
#include "graph_blockmodel_layers.hh"
#include "graph_blockmodel_multiflip_mcmc.hh"
#include "mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, BlockState, BLOCK_STATE_params)
template <class BaseState>
GEN_DISPATCH(layered_block_state, Layers<BaseState>::template LayeredBlockState,
LAYERED_BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object multiflip_mcmc_layered_sweep(python::object omcmc_state,
python::object olayered_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto* block_state)
{
typedef typename std::remove_pointer<decltype(block_state)>::type
state_t;
layered_block_state<state_t>::dispatch
(olayered_state,
[&](auto& ls)
{
typedef typename std::remove_reference<decltype(ls)>::type
layered_state_t;
mcmc_block_state<layered_state_t>::make_dispatch
(omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
});
},
false);
};
block_state::dispatch(dispatch);
return ret;
}
void export_layered_blockmodel_multiflip_mcmc()
{
using namespace boost::python;
def("multiflip_mcmc_layered_sweep", &multiflip_mcmc_layered_sweep);
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_blockmodel_overlap_util.hh"
#include "graph_blockmodel_overlap.hh"
#define BASE_STATE_params OVERLAP_BLOCK_STATE_params ((eweight,,,0))
#include "graph_blockmodel_layers.hh"
#include "graph_blockmodel_multiflip_mcmc.hh"
#include "mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(overlap_block_state, OverlapBlockState, OVERLAP_BLOCK_STATE_params)
template <class BaseState>
GEN_DISPATCH(layered_block_state, Layers<BaseState>::template LayeredBlockState,
LAYERED_BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object multiflip_mcmc_layered_overlap_sweep(python::object omcmc_state,
python::object olayered_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto* block_state)
{
typedef typename std::remove_pointer<decltype(block_state)>::type
state_t;
layered_block_state<state_t>::dispatch
(olayered_state,
[&](auto& ls)
{
typedef typename std::remove_reference<decltype(ls)>::type
layered_state_t;
mcmc_block_state<layered_state_t>::make_dispatch
(omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
});
},
false);
};
overlap_block_state::dispatch(dispatch);
return ret;
}
void export_layered_overlap_blockmodel_multiflip_mcmc()
{
using namespace boost::python;
def("multiflip_mcmc_layered_overlap_sweep",
&multiflip_mcmc_layered_overlap_sweep);
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_blockmodel_util.hh"
#include "graph_blockmodel.hh"
#include "graph_blockmodel_multiflip_mcmc.hh"
#include "mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, BlockState, BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object do_multiflip_mcmc_sweep(python::object omcmc_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto& block_state)
{
typedef typename std::remove_reference<decltype(block_state)>::type
state_t;
mcmc_block_state<state_t>::make_dispatch
(omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
});
};
block_state::dispatch(oblock_state, dispatch);
return ret;
}
void export_blockmodel_multiflip_mcmc()
{
using namespace boost::python;
def("multiflip_mcmc_sweep", &do_multiflip_mcmc_sweep);
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef GRAPH_BLOCKMODEL_MULTIFLIP_MCMC_HH
#define GRAPH_BLOCKMODEL_MULTIFLIP_MCMC_HH
#include "config.h"
#include <vector>
#include <algorithm>
#include "graph_tool.hh"
#include "graph_state.hh"
#include "graph_blockmodel_util.hh"
#include <boost/mpl/vector.hpp>
namespace graph_tool
{
using namespace boost;
using namespace std;
#define MCMC_BLOCK_STATE_params(State) \
((__class__,&, mpl::vector<python::object>, 1)) \
((state, &, State&, 0)) \
((E,, size_t, 0)) \
((beta,, double, 0)) \
((c,, double, 0)) \
((a,, double, 0)) \
((entropy_args,, entropy_args_t, 0)) \
((mproposals, & ,std::vector<size_t>&, 0)) \
((maccept, & ,std::vector<size_t>&, 0)) \
((allow_vacate,, bool, 0)) \
((sequential,, bool, 0)) \
((verbose,, bool, 0)) \
((niter,, size_t, 0))
template <class State>
struct MCMC
{
GEN_STATE_BASE(MCMCBlockStateBase, MCMC_BLOCK_STATE_params(State))
template <class... Ts>
class MCMCBlockState
: public MCMCBlockStateBase<Ts...>
{
public:
GET_PARAMS_USING(MCMCBlockStateBase<Ts...>,
MCMC_BLOCK_STATE_params(State))
GET_PARAMS_TYPEDEF(Ts, MCMC_BLOCK_STATE_params(State))
template <class... ATs,
typename std::enable_if_t<sizeof...(ATs) ==
sizeof...(Ts)>* = nullptr>
MCMCBlockState(ATs&&... as)
: MCMCBlockStateBase<Ts...>(as...),
_g(_state._g),
_groups(num_vertices(_state._bg)),
_vpos(get(vertex_index_t(), _state._g),
num_vertices(_state._g)),
_mprobs(num_vertices(_state._g) + 1)
{
_state.init_mcmc(_c,
(_entropy_args.partition_dl ||
_entropy_args.degree_dl ||
_entropy_args.edges_dl));
for (auto v : vertices_range(_state._g))
{
if (_state._vweight[v] > 0)
add_element(_groups[_state._b[v]], _vpos, v);
}
for (auto r : vertices_range(_state._bg))
_vlist.push_back(r);
for (size_t m = 1; m < _mprobs.size(); ++m)
_mprobs[m] = _mprobs[m-1] + std::pow(m, -_a);
}
typename state_t::g_t& _g;
std::vector<std::vector<size_t>> _groups;
typename vprop_map_t<size_t>::type::unchecked_t _vpos;
std::vector<size_t> _vlist;
std::vector<size_t> _vs;
std::vector<double> _mprobs;
size_t node_state(size_t r)
{
return r;
}
size_t node_weight(size_t r)
{
return _groups[r].size();
}
template <class RNG>
size_t sample_m(size_t n, RNG& rng)
{
double M = _mprobs[n];
std::uniform_real_distribution<> u_sample(0, M);
auto u = u_sample(rng);
auto iter = std::lower_bound(_mprobs.begin(),
_mprobs.begin() + n + 1, u);
return iter - _mprobs.begin();
}
double log_pm(size_t m, size_t n)
{
return - _a * log(m) - log(_mprobs[n]);
}
template <class RNG>
size_t move_proposal(size_t r, RNG& rng)
{
if (!_allow_vacate && _groups[r].size() == 1)
return r;
size_t m = sample_m(_groups[r].size(), rng);
assert(m <= _groups[r].size());
_vs.clear();
while (_vs.size() < m)
{
size_t v = uniform_sample(_groups[r], rng);
_vs.push_back(v);
remove_element(_groups[r], _vpos, v);
}
for (auto v : _vs)
add_element(_groups[r], _vpos, v);
size_t v = uniform_sample(_vs, rng);
size_t s = _state.sample_block(v, _c, rng);
if (!_state.allow_move(r, s))
return r;
if (_groups[s].size() == 0)
_state._bclabel[s] = _state._bclabel[r];
_mproposals[m]++;
return s;
}
std::tuple<double, double>
virtual_move_dS(size_t r, size_t nr)
{
double dS = 0, a = 0;
size_t m = _vs.size();
a -= log_pm(m, _groups[r].size());
a += log_pm(m, _groups[nr].size() + m);
a -= -lbinom(_groups[r].size(), m);
a += -lbinom(_groups[nr].size() + m, m);
if (m == 1)
{
auto v = _vs.front();
dS = _state.virtual_move(v, r, nr, _entropy_args);
if (!std::isinf(_c))
{
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);
}
}
else
{
_state._egroups_enabled = false;
double pf = 0, pb = 0;
if (!std::isinf(_c))
{
for (auto v : _vs)
pf += _state.get_move_prob(v, r, nr, _c, false);
pf /= m;
}
for (auto v : _vs)
{
dS += _state.virtual_move(v, r, nr, _entropy_args);
_state.move_vertex(v, nr);
}
if (!std::isinf(_c))
{
for (auto v : _vs)
pb += _state.get_move_prob(v, nr, r, _c, false);
pb /= m;
a += log(pb) - log(pf);
}
for (auto v : _vs)
_state.move_vertex(v, r);
_state._egroups_enabled = true;
}
return {dS, a};
}
void perform_move(size_t r, size_t nr)
{
for (auto v : _vs)
{
_state.move_vertex(v, nr);
remove_element(_groups[r], _vpos, v);
add_element(_groups[nr], _vpos, v);
}
_maccept[_vs.size()]++;
}
};
};
} // graph_tool namespace
#endif //GRAPH_BLOCKMODEL_MCMC_HH
......@@ -87,6 +87,7 @@ public:
_c_brec(_brec.get_checked()),
_c_bdrec(_bdrec.get_checked()),
_emat(_bg, rng),
_egroups_enabled(true),
_overlap_stats(_g, _b, _half_edges, _node_index, num_vertices(_bg))
{
for (auto r : vertices_range(_bg))
......@@ -101,6 +102,7 @@ public:
_c_brec(_brec.get_checked()),
_c_bdrec(_bdrec.get_checked()),
_emat(other._emat),
_egroups_enabled(other._egroups_enabled),
_overlap_stats(other._overlap_stats)
{
if (other.is_partition_stats_enabled())
......@@ -152,7 +154,7 @@ public:
_overlap_stats.remove_half_edge(v, r, _b, _g);
_wr[r] = _overlap_stats.get_block_size(r);
if (!_egroups.empty())
if (!_egroups.empty() && _egroups_enabled)
_egroups.remove_vertex(v, _b, _g);
}
......@@ -241,7 +243,7 @@ public:
_b[v] = r;
if (!_egroups.empty())
if (!_egroups.empty() && _egroups_enabled)
_egroups.add_vertex(v, _b, _eweight, _g);
}
......@@ -1012,6 +1014,7 @@ public:
emat_t _emat;
EGroups<g_t, mpl::false_> _egroups;
bool _egroups_enabled;
overlap_stats_t _overlap_stats;
std::vector<overlap_partition_stats_t> _partition_stats;
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_blockmodel_overlap_util.hh"
#include "graph_blockmodel_overlap.hh"
#include "graph_blockmodel_multiflip_mcmc.hh"
#include "mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(overlap_block_state, OverlapBlockState, OVERLAP_BLOCK_STATE_params)
template<class State>
using OverlapMCMC = MCMC<State>;
template <class State>
GEN_DISPATCH(mcmc_overlap_block_state,
OverlapMCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object overlap_multiflip_mcmc_sweep(python::object omcmc_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto& block_state)
{
typedef typename std::remove_reference<decltype(block_state)>::type
state_t;
mcmc_overlap_block_state<state_t>::make_dispatch
(omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
return ret;
}
void export_overlap_blockmodel_multiflip_mcmc()
{
using namespace boost::python;
def("overlap_multiflip_mcmc_sweep", &overlap_multiflip_mcmc_sweep);
}
......@@ -75,12 +75,14 @@ void vector_rmap(boost::python::object ovals, boost::python::object omap)
extern void export_blockmodel_state();
extern void export_blockmodel_mcmc();
extern void export_blockmodel_multicanonical();
extern void export_blockmodel_multiflip_mcmc();
extern void export_blockmodel_merge();
extern void export_blockmodel_gibbs();
extern void export_overlap_blockmodel_state();
extern void export_overlap_blockmodel_mcmc();
extern void export_overlap_blockmodel_mcmc_bundled();
extern void export_overlap_blockmodel_multicanonical();
extern void export_overlap_blockmodel_multiflip_mcmc();
extern void export_overlap_blockmodel_gibbs();
extern void export_overlap_blockmodel_vacate();
extern void export_layered_blockmodel_state();
......@@ -88,11 +90,13 @@ extern void export_layered_blockmodel_mcmc();