Commit e7446106 authored by Tiago Peixoto's avatar Tiago Peixoto

Implement reconstruction from dynamics

parent 6ea36ec4
......@@ -127,7 +127,7 @@ public:
_n_items--;
}
void clear(bool shrink)
void clear(bool shrink=false)
{
_items.clear();
_ipos.clear();
......
......@@ -55,6 +55,19 @@ libgraph_tool_inference_la_SOURCES = \
layers/graph_blockmodel_layers_overlap_multicanonical_multiflip.cc \
layers/graph_blockmodel_layers_overlap_multiflip_mcmc.cc \
layers/graph_blockmodel_layers_overlap_vacate.cc \
uncertain/graph_blockmodel_dynamics_epidemics.cc \
uncertain/graph_blockmodel_dynamics_epidemics_mcmc.cc \
uncertain/graph_blockmodel_dynamics_epidemics_mcmc_r.cc \
uncertain/graph_blockmodel_dynamics_cising_glauber.cc \
uncertain/graph_blockmodel_dynamics_cising_glauber_mcmc.cc \
uncertain/graph_blockmodel_dynamics_ising_glauber.cc \
uncertain/graph_blockmodel_dynamics_ising_glauber_mcmc.cc \
uncertain/graph_blockmodel_dynamics_pseudo_cising.cc \
uncertain/graph_blockmodel_dynamics_pseudo_cising_mcmc.cc \
uncertain/graph_blockmodel_dynamics_pseudo_cising_mcmc_h.cc \
uncertain/graph_blockmodel_dynamics_pseudo_ising.cc \
uncertain/graph_blockmodel_dynamics_pseudo_ising_mcmc.cc \
uncertain/graph_blockmodel_dynamics_pseudo_ising_mcmc_h.cc \
uncertain/graph_blockmodel_measured.cc \
uncertain/graph_blockmodel_measured_mcmc.cc \
uncertain/graph_blockmodel_uncertain.cc \
......@@ -89,6 +102,11 @@ libgraph_tool_inference_la_include_HEADERS = \
overlap/graph_blockmodel_overlap_partition.hh \
layers/graph_blockmodel_layers.hh \
layers/graph_blockmodel_layers_util.hh \
uncertain/graph_blockmodel_dynamics.hh \
uncertain/graph_blockmodel_dynamics_discrete.hh \
uncertain/graph_blockmodel_dynamics_continuous.hh \
uncertain/graph_blockmodel_dynamics_epidemics_mcmc_h.hh \
uncertain/graph_blockmodel_dynamics_pseudo_ising_mcmc_h.hh \
uncertain/graph_blockmodel_measured.hh \
uncertain/graph_blockmodel_uncertain.hh \
uncertain/graph_blockmodel_uncertain_marginal.hh \
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2018 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_MULTILEVEL_MCMC_HH
#define GRAPH_BLOCKMODEL_MULTILEVEL_MCMC_HH
#include "config.h"
#include <vector>
#include <algorithm>
#include "graph_tool.hh"
#include "../support/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)) \
((beta,, double, 0)) \
((c,, double, 0)) \
((a,, double, 0)) \
((entropy_args,, entropy_args_t, 0)) \
((allow_vacate,, 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)),
_rpos(get(vertex_index_t(), _state._bg),
num_vertices(_state._bg)),
_clabel(get(vertex_index_t(), _state._g),
num_vertices(_state._g))
{
_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)
{
auto r = _state._b[v];
if (_groups[r].empty())
add_element(_rlist, _rpos, r);
add_element(_groups[r], _vpos, v);
}
}
}
typename state_t::g_t& _g;
std::vector<std::vector<size_t>> _groups;
typename vprop_map_t<size_t>::type::unchecked_t _vpos;
typename vprop_map_t<size_t>::type::unchecked_t _rpos;
std::vector<size_t> _rlist;
typename vprop_map_t<size_t>::type::unchecked_t _blabel;
std::vector<gt_hash_map<size_t, size_t>> _mrpc;
size_t get_mrpc(size_t r, size_t c)
{
auto& mrp = _mrpc[r];
auto iter = mrp.find(c);
if (iter != mrp.end())
return iter->second;
return 0;
}
std::vector<size_t> seed_rs;
typename vprop_map_t<size_t>::type::unchecked_t _rtree;
std::vector<std::pair<size_t, size_t>> _merges;
template <class RNG>
size_t sample_merge(size_t r, RNG& rng)
{
if (_state._mrp[r] + _state._mrm[r] == 0)
return uniform_sample(_rclist[_blabel[r]], rng);
const auto& e = _egroups.sample_edge(r, rng);
auto t = _state._b[target(e, _g)];
if (t == r)
t = _b[source(e, _g)];
double p_rand = 0;
if (c > 0)
{
size_t B = _rlist.size();
auto etr = _state.get_mrs(t, r);
if (graph_tool::is_directed(_g))
etr += _state.get_mrs(r, t);
p_rand = c * B / double(get_mrpc(t, _blabel[r]) - etr + c * B);
}
size_t s;
std::uniform_real_distribution<> rdist;
if (prand > 0 && rdist(rng) < p_rand)
{
s = uniform_sample(_rclist[_blabel[r]], rng);
}
else
{
do
{
const auto& e = _egroups.sample_edge(t, rng);
auto s = _state._b[target(e, _g)];
if (s == t)
s = _b[source(e, _g)];
}
while (s == r || _blabel[s] != _blabel[r]);
}
return s;
}
template <class RNG>
void build_tree(RNG& rng)
{
size_t m = sample_m(_groups[r].size(), rng);
if (!_allow_vacate && _groups[r].size() == m)
return null_group;
assert(m <= _groups[r].size());
_nproposals += m;
_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, _d, rng);
if (s >= _groups.size())
{
_groups.resize(s + 1);
_rpos.resize(s + 1);
}
if (!_state.allow_move(r, s) || s == r)
return null_group;
if (!_groups[s].empty() || _groups[r].size() > m)
_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);
int B = _vlist.size();
a -= -log(B);
if (_groups[r].size() == m)
B--;
if (_groups[nr].empty())
B++;
a += -log(B);
if (m == 1)
{
auto v = _vs.front();
dS = _state.virtual_move(v, r, nr, _entropy_args);
double pf = _state.get_move_prob(v, r, nr, _c, _d, false);
double pb = _state.get_move_prob(v, nr, r, _c, _d, true);
a += log(pb) - log(pf);
}
else
{
_state._egroups_enabled = false;
double pf = 0, pb = 0;
for (auto v : _vs)
pf += _state.get_move_prob(v, r, nr, _c, _d, false);
pf /= m;
for (auto v : _vs)
{
dS += _state.virtual_move(v, r, nr, _entropy_args);
_state.move_vertex(v, nr);
}
for (auto v : _vs)
pb += _state.get_move_prob(v, nr, r, _c, _d, false);
pb /= m;
a += log(pb) - log(pf);
for (auto v : _vs)
_state.move_vertex(v, r);
_state._egroups_enabled = true;
}
return std::make_tuple(dS, a);
}
void perform_move(size_t r, size_t nr)
{
if (!_groups[nr].empty() || _groups[r].size() > _vs.size())
_maccept[_vs.size()]++;
if (_state._wr[nr] == 0)
add_element(_vlist, _rpos, nr);
for (auto v : _vs)
{
_state.move_vertex(v, nr);
remove_element(_groups[r], _vpos, v);
add_element(_groups[nr], _vpos, v);
}
if (_state._wr[r] == 0)
remove_element(_vlist, _rpos, r);
}
bool is_deterministic()
{
return false;
}
bool is_sequential()
{
return false;
}
auto& get_vlist()
{
return _vlist;
}
double get_beta()
{
return _beta;
}
size_t get_niter()
{
if (_nproposals < _niter * _N)
return std::numeric_limits<size_t>::max();
return 0;
}
void step(size_t, size_t)
{
}
};
};
} // graph_tool namespace
#endif //GRAPH_BLOCKMODEL_MCMC_HH
......@@ -108,8 +108,21 @@ extern void export_uncertain_state();
extern void export_uncertain_mcmc();
extern void export_measured_state();
extern void export_measured_mcmc();
extern void export_epidemics_state();
extern void export_epidemics_mcmc();
extern void export_epidemics_mcmc_r();
extern void export_cising_glauber_state();
extern void export_cising_glauber_mcmc();
extern void export_ising_glauber_state();
extern void export_ising_glauber_mcmc();
extern void export_marginals();
extern void export_modularity();
extern void export_pseudo_cising_state();
extern void export_pseudo_cising_mcmc();
extern void export_pseudo_cising_mcmc_h();
extern void export_pseudo_ising_state();
extern void export_pseudo_ising_mcmc();
extern void export_pseudo_ising_mcmc_h();
BOOST_PYTHON_MODULE(libgraph_tool_inference)
{
......@@ -154,8 +167,21 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
export_uncertain_mcmc();
export_measured_state();
export_measured_mcmc();
export_epidemics_state();
export_epidemics_mcmc();
export_epidemics_mcmc_r();
export_cising_glauber_state();
export_cising_glauber_mcmc();
export_ising_glauber_state();
export_ising_glauber_mcmc();
export_marginals();
export_modularity();
export_pseudo_cising_state();
export_pseudo_cising_mcmc();
export_pseudo_cising_mcmc_h();
export_pseudo_ising_state();
export_pseudo_ising_mcmc();
export_pseudo_ising_mcmc_h();
def("vector_map", vector_map<int32_t>);
def("vector_map64", vector_map<int64_t>);
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2018 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_DYNAMICS_HH
#define GRAPH_BLOCKMODEL_DYNAMICS_HH
#include "config.h"
#include <vector>
#include "../support/graph_state.hh"
#include "../blockmodel/graph_blockmodel_util.hh"
#include "graph_blockmodel_uncertain_util.hh"
namespace graph_tool
{
using namespace boost;
using namespace std;
typedef typename eprop_map_t<double>::type xmap_t;
#define DYNAMICS_STATE_params \
((__class__,&, mpl::vector<python::object>, 1)) \
((g, &, all_graph_views, 1)) \
((params,, python::dict, 0)) \
((ot,, python::list, 0)) \
((os,, python::list, 0)) \
((x,, xmap_t, 0)) \
((aE,, double, 0)) \
((E_prior,, bool, 0)) \
((self_loops,, bool, 0))
template <class Type>
std::vector<typename Type::unchecked_t> from_list(python::object list)
{
vector<typename Type::unchecked_t> v;
for (int i = 0; i < python::len(list); ++i)
{
boost::any a = python::extract<boost::any>(list[i])();
v.push_back(boost::any_cast<Type>(a).get_unchecked());
}
return v;
};
template <class BlockState, class DState>
struct Dynamics
{
GEN_STATE_BASE(DynamicsStateBase, DYNAMICS_STATE_params)
template <class... Ts>
class DynamicsState
: public DynamicsStateBase<Ts...>
{
public:
GET_PARAMS_USING(DynamicsStateBase<Ts...>,
DYNAMICS_STATE_params)
GET_PARAMS_TYPEDEF(Ts, DYNAMICS_STATE_params)
typedef typename DState::tmap_t tmap_t;
typedef typename DState::smap_t smap_t;
template <class... ATs,
typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
DynamicsState(BlockState& block_state, ATs&&... args)
: DynamicsStateBase<Ts...>(std::forward<ATs>(args)...),
_block_state(block_state),
_t(from_list<tmap_t>(_ot)),
_s(from_list<smap_t>(_os)),
_dstate(*this, _params),
_xc(_x.get_checked())
{
_u_edges.resize(num_vertices(_u));
for (auto e : edges_range(_u))
{
get_u_edge<true>(source(e, _u), target(e, _u)) = e;
_E += _eweight[e];
}
_block_state.enable_partition_stats();
}
DynamicsState(const DynamicsState& other)
: DynamicsStateBase<Ts...>(static_cast<const DynamicsStateBase<Ts...>&>(other)),
_block_state(other._block_state),
_t(other._t),
_s(other._s),
_u_edges(other._u_edges),
_pe(other._pe),
_E(other._E),
_dstate(*this, _params),
_xc(_x.get_checked())
{
_block_state.enable_partition_stats();
}
typedef BlockState block_state_t;
BlockState& _block_state;
std::vector<typename tmap_t::unchecked_t> _t;
std::vector<typename smap_t::unchecked_t> _s;
typename BlockState::g_t& _u = _block_state._g;
typename BlockState::eweight_t& _eweight = _block_state._eweight;
GraphInterface::edge_t _null_edge;
std::vector<double> _recs;
std::vector<gt_hash_map<size_t, GraphInterface::edge_t>> _u_edges;
double _pe = log(_aE);
size_t _E = 0;
DState _dstate;
xmap_t _xc;
template <bool insert, class Graph, class Elist>
auto& _get_edge(size_t u, size_t v, Graph& g, Elist& edges)
{
if (!graph_tool::is_directed(g) && u > v)
std::swap(u, v);
auto& qe = edges[u];
if (insert)
return qe[v];
auto iter = qe.find(v);
if (iter != qe.end())
return iter->second;
return _null_edge;
}
template <bool insert=false>
auto& get_u_edge(size_t u, size_t v)
{
return _get_edge<insert>(u, v, _u, _u_edges);
}
double get_node_prob(size_t u)
{
return _dstate.get_node_prob(u);
}
double entropy(bool latent_edges, bool density)
{
double S = 0;
if (latent_edges)
{
for (auto v : vertices_range(_u))
S += _dstate.get_node_prob(v);
}
if (density && _E_prior)
S += _E * _pe - lgamma_fast(_E + 1) - exp(_pe);
return -S;
}
double remove_edge_dS(size_t u, size_t v, uentropy_args_t ea)
{
auto& e = get_u_edge(u, v);
auto x = _xc[e];
double dS = _block_state.template modify_edge_dS<false>(source(e, _u),
target(e, _u),
e, _recs, ea);
_xc[e] = x;
if (ea.density && _E_prior)
{
dS += _pe;
dS += lgamma_fast(_E) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if (_eweight[e] == 1 && (_self_loops || u != v))
{
dS += _dstate.template get_edge_dS<false>(u, v, _xc[e]);
if (u != v && !graph_tool::is_directed(_u))
dS += _dstate.template get_edge_dS<false>(v, u, _xc[e]);
}
}
return dS;
}
double add_edge_dS(size_t u, size_t v, double x, uentropy_args_t ea)
{
auto& e = get_u_edge(u, v);
double dS = _block_state.template modify_edge_dS<true>(u, v, e,
_recs, ea);
if (ea.density && _E_prior)
{
dS -= _pe;
dS += lgamma_fast(_E + 2) - lgamma_fast(_E + 1);
}
if (ea.latent_edges)
{
if ((e == _null_edge || _eweight[e] == 0) && (_self_loops || u != v))
{
dS += _dstate.template get_edge_dS<true>(u, v, x);
if (u != v && !graph_tool::is_directed(_u))
dS += _dstate.template get_edge_dS<true>(v, u, x);
}
}
return dS;
}
double update_edge_dS(size_t u, size_t v, double dx, uentropy_args_t ea)
{
double dS = 0;
if (ea.latent_edges)
{
if (_self_loops || u != v)
{
dS += _dstate.template get_edge_dS<true>(u, v, dx);
if (u != v && !graph_tool::is_directed(_u))
dS += _dstate.template get_edge_dS<true>(v, u, dx);
}
}
return dS;
}
void remove_edge(size_t u, size_t v)
{
auto& e = get_u_edge(u, v);
auto x = _xc[e];
_block_state.template modify_edge<false>(u, v, e,
_recs);
if ((e == _null_edge || _eweight[e] == 0) && (_self_loops || u != v))
{
_dstate.template update_edge<false>(u, v, x);
if (u != v && !graph_tool::is_directed(_u))
_dstate.template update_edge<false>(v, u, x);
}
_E--;
}
void add_edge(size_t u, size_t v, double x)
{
auto& e = get_u_edge<true>(u, v);
_block_state.template modify_edge<true>(u, v, e,
_recs);
if (_eweight[e] == 1 && (_self_loops || u != v))
{
_xc[e] = x;
_dstate.template update_edge<true>(u, v, x);
if (u != v && !graph_tool::is_directed(_u))
_dstate.template update_edge<true>(v, u, x);
}
_E++;
}
void update_edge(size_t u, size_t v, double dx)
{
if (_self_loops || u != v)
{
auto& e = get_u_edge(u, v);
_xc[e] += dx;
_dstate.template update_edge<true>(u, v, dx);
if (u != v && !graph_tool::is_directed(_u))
_dstate.template update_edge<true>(v, u, dx);
}
}
void set_params(python::dict params)
{
_dstate.set_params(params);
}
};
};
} // graph_tool namespace
#endif //GRAPH_BLOCKMODEL_DYNAMICS_HH