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

Implement Ranked SBM

parent e24b5387
......@@ -425,6 +425,11 @@ libgraph_tool_inference_la_SOURCES = \
src/graph/inference/layers/graph_blockmodel_layers_overlap_multiflip_mcmc.cc \
src/graph/inference/layers/graph_blockmodel_layers_overlap_multilevel_mcmc.cc \
src/graph/inference/layers/graph_blockmodel_layers_overlap_vacate.cc \
src/graph/inference/ranked/graph_ranked.cc \
src/graph/inference/ranked/graph_ranked_mcmc.cc \
src/graph/inference/ranked/graph_ranked_gibbs.cc \
src/graph/inference/ranked/graph_ranked_multiflip_mcmc.cc \
src/graph/inference/ranked/graph_ranked_multilevel_mcmc.cc \
src/graph/inference/uncertain/graph_blockmodel_dynamics_epidemics.cc \
src/graph/inference/uncertain/graph_blockmodel_dynamics_epidemics_mcmc.cc \
src/graph/inference/uncertain/graph_blockmodel_dynamics_epidemics_mcmc_r.cc \
......@@ -501,6 +506,14 @@ graph_inference_planted_partitiondir = $(MOD_DIR)/include/inference/planted_part
dist_graph_inference_planted_partition_HEADERS = \
src/graph/inference/planted_partition/graph_planted_partition.hh
graph_inference_rankeddir = $(MOD_DIR)/include/inference/ranked
dist_graph_inference_ranked_HEADERS = \
src/graph/inference/ranked/graph_ranked.hh \
src/graph/inference/ranked/graph_ranked_mcmc.hh \
src/graph/inference/ranked/graph_ranked_gibbs.hh \
src/graph/inference/ranked/graph_ranked_multiflip_mcmc.hh \
src/graph/inference/ranked/graph_ranked_multilevel_mcmc.hh
graph_inference_overlapdir = $(MOD_DIR)/include/inference/overlap
dist_graph_inference_overlap_HEADERS = \
src/graph/inference/overlap/graph_blockmodel_overlap.hh \
......@@ -787,6 +800,7 @@ graph_tool_inference_PYTHON = \
src/graph_tool/inference/partition_centroid.py \
src/graph_tool/inference/partition_modes.py \
src/graph_tool/inference/planted_partition.py \
src/graph_tool/inference/ranked.py \
src/graph_tool/inference/latent_multigraph.py \
src/graph_tool/inference/latent_layers.py \
src/graph_tool/inference/util.py
......
......@@ -17,6 +17,7 @@ reader is referred to [peixoto-bayesian-2019]_.
.. include:: _edge_weights.rst
.. include:: _layers.rst
.. include:: _assortative.rst
.. include:: _ranked.rst
.. include:: _reconstruction.rst
.. include:: _prediction.rst
......@@ -86,6 +87,9 @@ References
inference of assortative community structures", Phys. Rev. Research 2
043271 (2020), :doi:`10.1103/PhysRevResearch.2.043271`, :arxiv:`2006.14493`
.. [peixoto-ordered-2022] Tiago P. Peixoto, "Ordered community detection
in directed networks", :arxiv:`2203.16460`
.. [peixoto-reconstructing-2018] Tiago P. Peixoto, "Reconstructing
networks with unknown and heterogeneous errors", Phys. Rev. X 8
041011 (2018). :doi:`10.1103/PhysRevX.8.041011`, :arxiv:`1806.07956`
......
......@@ -144,6 +144,11 @@ extern void export_pp_state();
extern void export_pp_mcmc();
extern void export_pp_multiflip_mcmc();
extern void export_pp_multilevel_mcmc();
extern void export_ranked_gibbs();
extern void export_ranked_state();
extern void export_ranked_mcmc();
extern void export_ranked_multiflip_mcmc();
extern void export_ranked_multilevel_mcmc();
extern void export_modularity_state();
extern void export_modularity_gibbs();
extern void export_modularity_mcmc();
......@@ -239,6 +244,11 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
export_pp_mcmc();
export_pp_multiflip_mcmc();
export_pp_multilevel_mcmc();
export_ranked_state();
export_ranked_gibbs();
export_ranked_mcmc();
export_ranked_multiflip_mcmc();
export_ranked_multilevel_mcmc();
export_modularity_state();
export_modularity_gibbs();
export_modularity_mcmc();
......
......@@ -929,9 +929,7 @@ struct MergeSplit: public State
{
auto r = uniform_sample(_rlist, rng);
auto s = r;
while (get_wr(s) > 0)
s = sample_move(r, rng);
auto s = State::sample_new_group(*_groups[r].begin(), rng);
if (!allow_merge(r, s))
{
......@@ -941,9 +939,6 @@ struct MergeSplit: public State
State::relax_update(true);
if (!std::isinf(_beta))
pf += merge_prob(r, s); //FIXME
auto& vrs = _groups[r];
_vs.clear();
_vs.insert(_vs.begin(), vrs.begin(), vrs.end());
......@@ -955,9 +950,6 @@ struct MergeSplit: public State
for (const auto& v : _vs)
_bnext[v] = State::get_group(v);
if (!std::isinf(_beta))
pb += merge_prob(s, r); //FIXME
pop_b();
State::relax_update(false);
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser 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 "../support/graph_state.hh"
#define GRAPH_RANGE never_filtered_never_reversed
#include "../blockmodel/graph_blockmodel.hh"
#include "graph_ranked.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, BlockState, BLOCK_STATE_params)
template <class BState>
GEN_DISPATCH(ranked_state, OState<BState>::template RankedState, RANKED_STATE_params)
python::object make_ranked_state(boost::python::object oustate,
boost::python::object ostate)
{
python::object state;
auto dispatch_d = [&](auto& ustate)
{
typedef typename std::remove_reference<decltype(ustate)>::type
ustate_t;
ranked_state<ustate_t>
::make_dispatch(ostate,
[&](auto& s){state = python::object(s);},
ustate);
};
block_state::dispatch(oustate, dispatch_d);
return state;
}
void export_ranked_state()
{
using namespace boost::python;
def("make_ranked_state", &make_ranked_state);
block_state::dispatch
([&](auto* ds)
{
typedef typename std::remove_reference<decltype(*ds)>::type dstate_t;
ranked_state<dstate_t>::dispatch
([&](auto* s)
{
typedef typename std::remove_reference<decltype(*s)>::type state_t;
void (state_t::*move_vertex)(size_t, size_t) =
&state_t::move_vertex;
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t& ea) =
&state_t::virtual_move;
class_<state_t, bases<>, std::shared_ptr<state_t>>
c(name_demangle(typeid(state_t).name()).c_str(),
no_init);
c.def("move_vertex", move_vertex)
.def("virtual_move", virtual_move)
.def("entropy", &state_t::entropy)
.def("couple_state", &state_t::couple_state)
.def("decouple_state", &state_t::decouple_state)
.def("get_Es",
+[](state_t& state)
{
return python::make_tuple(state._E[0],
state._E[1],
state._E[2]);
});
});
});
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef GRAPH_RANKED_HH
#define GRAPH_RANKED_HH
#include "config.h"
#include <vector>
#include "../blockmodel/graph_blockmodel_util.hh"
#include "../support/graph_state.hh"
namespace graph_tool
{
using namespace boost;
using namespace std;
typedef vprop_map_t<int32_t>::type vmap_t;
typedef eprop_map_t<int32_t>::type emap_t;
typedef vprop_map_t<double>::type umap_t;
#define RANKED_STATE_params \
((__class__,&, mpl::vector<python::object>, 1)) \
((u,, umap_t, 0))
GEN_STATE_BASE(RankedStateBase, RANKED_STATE_params)
template <class BState>
class OState
{
public:
template <class... Ts>
class RankedState
: public RankedStateBase<Ts...>
{
public:
GET_PARAMS_USING(RankedStateBase<Ts...>, RANKED_STATE_params)
GET_PARAMS_TYPEDEF(Ts, RANKED_STATE_params)
typedef BState bstate_t;
typedef typename bstate_t::_entropy_args_t _entropy_args_t;
template <class... ATs,
typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
RankedState(BState& ustate, ATs&&... args)
: RankedStateBase<Ts...>(std::forward<ATs>(args)...),
_ustate(ustate),
_g(ustate._g),
_b(ustate._b),
_eweight(ustate._eweight),
_u_c(_u.get_checked())
{
for (auto e : edges_range(_g))
{
auto u = source(e, _g);
auto v = target(e, _g);
_E[stream_dir(_b[u], _b[v])] += _eweight[e];
}
}
bstate_t& _ustate;
typename bstate_t::g_t& _g;
typename bstate_t::b_t& _b;
typename bstate_t::eweight_t& _eweight;
std::array<size_t, 3> _E = {0, 0, 0};
typename u_t::checked_t _u_c;
int stream_dir(size_t r, size_t s)
{
auto x = _u[r];
auto y = _u[s];
if (x < y)
return 0; // upstream
if (x > y)
return 2; // downstream
return 1;
}
void add_block(size_t n)
{
_ustate.add_block(n);
}
gt_hash_map<size_t, int> _delta;
double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t& ea)
{
entropy_args_t uea = ea;
uea.edges_dl = false;
double dS = _ustate.virtual_move(v, r, nr, uea);
std::array<int, 3> dE = {0, 0, 0};
dS -= get_edges_dl(dE, 0);
for (auto e : out_edges_range(v, _g))
{
auto s = _b[target(e, _g)];
auto w = _eweight[e];
dE[stream_dir(r, s)] -= w;
if (target(e, _g) == v)
s = nr;
dE[stream_dir(nr, s)] += w;
}
for (auto e : in_edges_range(v, _g))
{
auto s = _b[source(e, _g)];
auto w = _eweight[e];
dE[stream_dir(s, r)] -= w;
if (source(e, _g) == v)
s = nr;
dE[stream_dir(s, nr)] += w;
}
int dB = 0;
if (_ustate._wr[r] == 1)
dB--;
if (_ustate._wr[nr] == 0)
dB++;
dS += get_edges_dl(dE, dB);
_delta.clear();
size_t B = num_vertices(_ustate._bg) + 1;
entries_op(_ustate._m_entries, _ustate._emat,
[&](auto t, auto u, auto&, auto delta)
{
if (delta == 0 || t == u)
return;
_delta[t + B * u] = delta;
});
entries_op(_ustate._m_entries, _ustate._emat,
[&](auto t, auto u, auto& me, auto delta)
{
if (delta == 0 || t == u)
return;
size_t etu = _ustate._mrs[me];
size_t eut = get_beprop(u, t, _ustate._mrs,
_ustate._emat);
int delta_r = 0;
auto iter = _delta.find(u + B * t);
if (iter != _delta.end())
delta_r = iter->second;
if (delta_r != 0 && t > u)
return;
dS += lbinom_fast(etu + eut, etu);
dS -= lbinom_fast(etu + delta + eut + delta_r,
etu + delta);
});
return dS;
}
void move_vertex(size_t v, size_t nr)
{
auto r = _b[v];
for (auto e : out_edges_range(v, _g))
{
auto s = _b[target(e, _g)];
auto w = _eweight[e];
_E[stream_dir(r, s)] -= w;
if (target(e, _g) == v)
s = nr;
_E[stream_dir(nr, s)] += w;
}
for (auto e : in_edges_range(v, _g))
{
auto s = _b[source(e, _g)];
auto w = _eweight[e];
_E[stream_dir(s, r)] -= w;
if (source(e, _g) == v)
s = nr;
_E[stream_dir(s, nr)] += w;
}
_ustate.move_vertex(v, nr);
}
size_t get_empty_block(size_t v, bool force_add = false)
{
return _ustate.get_empty_block(v, force_add);
}
double sample_u(rng_t& rng)
{
uniform_real_distribution<double> usample;
return usample(rng);
}
size_t sample_block(size_t v, double c, double d, rng_t& rng)
{
auto s = _ustate.sample_block(v, c, d, rng);
if (_ustate._wr[s] == 0)
{
uniform_real_distribution<double> usample;
_u_c[s] = usample(rng);
}
return s;
}
size_t sample_block_local(size_t v, rng_t& rng)
{
return _ustate.sample_block_local(v, rng);
}
// Computes the move proposal probability
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse)
{
return _ustate.get_move_prob(v, r, s, c, d, reverse);
}
double get_edges_dl(const std::array<int,3>& dE, int dB)
{
size_t B = _ustate._candidate_groups.size() + dB;
double S = 0;
if (_ustate._coupled_state == nullptr)
S += lbinom_fast<false>((B * (B + 1)) / 2 + _ustate._E - 1,
_ustate._E);
std::array<size_t, 3> E;
for (int i = 0; i < 3; ++i)
E[i] = _E[i] + dE[i];
S += lgamma_fast(E[0] + E[2] + 2);
S -= lgamma_fast(E[0] + 1) + lgamma_fast(E[2] + 1);
return S;
}
double entropy(entropy_args_t ea)
{
double S = 0;
ea.edges_dl = false;
S += _ustate.entropy(ea);
S += get_edges_dl({0, 0, 0}, 0);
for (auto e : edges_range(_ustate._bg))
{
auto r = source(e, _ustate._bg);
auto s = target(e, _ustate._bg);
if (r >= s)
continue;
size_t ers = _ustate._mrs[e];
size_t esr = get_beprop(s, r, _ustate._mrs, _ustate._emat);
S -= lbinom_fast(ers + esr, ers);
}
return S;
}
template <class MCMCState>
void init_mcmc(MCMCState& state)
{
_ustate.init_mcmc(state);
}
size_t node_weight(size_t v)
{
return _ustate.node_weight(v);
}
bool is_last(size_t v)
{
return _ustate.is_last(v);
}
bool allow_move(size_t v, size_t r)
{
return _ustate.allow_move(v, r);
}
size_t virtual_remove_size(size_t v)
{
return _ustate.virtual_remove_size(v);
}
template <class V>
void push_state(V&) {}
void pop_state() {}
void store_next_state(size_t) {}
void clear_next_state() {}
void relax_update(bool relax)
{
_ustate.relax_update(relax);
}
void couple_state(BlockStateVirtualBase& us,
const entropy_args_t& ea)
{
_ustate.couple_state(us, ea);
}
void decouple_state()
{
_ustate.decouple_state();
}
};
};
} // graph_tool namespace
#endif //GRAPH_RANKED_HH
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser 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>
#define GRAPH_RANGE never_filtered_never_reversed
#include "../blockmodel/graph_blockmodel.hh"
#include "graph_ranked.hh"
#include "graph_ranked_gibbs.hh"
#include "../loops/gibbs_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, BlockState, BLOCK_STATE_params)
template <class BState>
GEN_DISPATCH(ranked_state, OState<BState>::template RankedState, RANKED_STATE_params)
template <class State>
GEN_DISPATCH(gibbs_block_state, Gibbs<State>::template GibbsBlockState,
GIBBS_BLOCK_STATE_params(State))
python::object ranked_gibbs_sweep(python::object ogibbs_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;
auto dispatch_d = [&](auto* dstate)
{
typedef typename std::remove_reference<decltype(*dstate)>::type
dstate_t;
auto dispatch = [&](auto& block_state)
{
typedef typename std::remove_reference<decltype(block_state)>::type
state_t;
gibbs_block_state<state_t>::make_dispatch
(ogibbs_state,
[&](auto& s)
{