Commit 100f01d4 authored by Tiago Peixoto's avatar Tiago Peixoto

Implement partition modes, centroid and manipulation functions

parent d5f2e2c9
Pipeline #617 passed with stage
in 825 minutes and 22 seconds
......@@ -64,6 +64,14 @@
:no-undoc-members:
:show-inheritance:
.. automodule:: graph_tool.inference.partition_modes
:no-undoc-members:
:show-inheritance:
.. automodule:: graph_tool.inference.partition_centroid
:no-undoc-members:
:show-inheritance:
.. automodule:: graph_tool.inference.blockmodel_em
:no-undoc-members:
:show-inheritance:
......
......@@ -31,6 +31,9 @@ libgraph_tool_inference_la_SOURCES = \
partition_centroid/graph_partition_centroid.cc \
partition_centroid/graph_partition_centroid_mcmc.cc \
partition_centroid/graph_partition_centroid_multiflip_mcmc.cc \
partition_centroid/graph_partition_centroid_rmi.cc \
partition_centroid/graph_partition_centroid_rmi_mcmc.cc \
partition_centroid/graph_partition_centroid_rmi_multiflip_mcmc.cc \
partition_modes/graph_partition_mode.cc \
partition_modes/graph_partition_mode_clustering.cc \
partition_modes/graph_partition_mode_clustering_mcmc.cc \
......@@ -104,6 +107,7 @@ nobase_libgraph_tool_inference_la_include_HEADERS = \
blockmodel/graph_blockmodel_util.hh \
blockmodel/graph_blockmodel_weights.hh \
partition_centroid/graph_partition_centroid.hh \
partition_centroid/graph_partition_centroid_rmi.hh \
partition_modes/graph_partition_mode.hh \
partition_modes/graph_partition_mode_clustering.hh \
overlap/graph_blockmodel_overlap.hh \
......
......@@ -508,10 +508,8 @@ public:
return S_a - S_b;
}
template <class Graph, class VWeight, class EWeight, class Degs>
void change_vertex(size_t v, size_t r, bool deg_corr, Graph& g,
VWeight& vweight, EWeight& eweight, Degs& degs,
int diff)
template <class VWeight>
void change_vertex(size_t v, size_t r, VWeight& vweight, int diff)
{
int vw = vweight[v];
int dv = vw * diff;
......@@ -526,23 +524,25 @@ public:
_N += dv;
assert(_total[r] >= 0);
}
if (deg_corr)
{
degs_op(v, vweight, eweight, degs, g,
[&](auto kin, auto kout, auto n)
{
int dk = diff * n;
auto& h = _hist[r];
auto deg = make_pair(kin, kout);
auto iter = h.insert({deg, 0}).first;
iter->second += dk;
if (iter->second == 0)
h.erase(iter);
_em[r] += dk * deg.first;
_ep[r] += dk * deg.second;
});
}
template <class Graph, class VWeight, class EWeight, class Degs>
void change_vertex_degs(size_t v, size_t r, Graph& g, VWeight& vweight,
EWeight& eweight, Degs& degs, int diff)
{
degs_op(v, vweight, eweight, degs, g,
[&](auto kin, auto kout, auto n)
{
int dk = diff * n;
auto& h = _hist[r];
auto deg = make_pair(kin, kout);
auto iter = h.insert({deg, 0}).first;
iter->second += dk;
if (iter->second == 0)
h.erase(iter);
_em[r] += dk * deg.first;
_ep[r] += dk * deg.second;
});
}
template <class Graph, class VWeight, class EWeight, class Degs>
......@@ -552,7 +552,9 @@ public:
if (r == null_group || vweight[v] == 0)
return;
r = get_r(r);
change_vertex(v, r, deg_corr, g, vweight, eweight, degs, -1);
change_vertex(v, r, vweight, -1);
if (deg_corr)
change_vertex_degs(v, r, g, vweight, eweight, degs, -1);
}
template <class Graph, class VWeight, class EWeight, class Degs>
......@@ -562,7 +564,9 @@ public:
if (nr == null_group || vweight[v] == 0)
return;
nr = get_r(nr);
change_vertex(v, nr, deg_corr, g, vweight, eweight, degs, 1);
change_vertex(v, nr, vweight, 1);
if (deg_corr)
change_vertex_degs(v, nr, g, vweight, eweight, degs, 1);
}
void change_k(size_t v, size_t r, bool deg_corr, int vweight,
......
......@@ -127,6 +127,9 @@ extern void export_pseudo_ising_mcmc_h();
extern void export_vi_center_state();
extern void export_vi_center_mcmc();
extern void export_vi_multiflip_mcmc();
extern void export_rmi_center_state();
extern void export_rmi_center_mcmc();
extern void export_rmi_multiflip_mcmc();
extern void export_partition_mode();
extern void export_mode_cluster_state();
extern void export_mode_cluster_mcmc();
......@@ -194,6 +197,9 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
export_vi_center_state();
export_vi_center_mcmc();
export_vi_multiflip_mcmc();
export_rmi_center_state();
export_rmi_center_mcmc();
export_rmi_multiflip_mcmc();
export_partition_mode();
export_mode_cluster_state();
export_mode_cluster_mcmc();
......
......@@ -129,6 +129,7 @@ public:
_wr[r]--;
_wr[nr]++;
//#pragma omp parallel for schedule(runtime)
for (size_t i = 0; i < _mrs.size(); ++i)
{
auto& mrsi = _mrs[i];
......@@ -173,6 +174,11 @@ public:
double Sb = 0;
double Sa = 0;
Sb += (xlogx_fast(_wr[r]) + xlogx_fast(_wr[nr])) * _mrs.size();
Sa += (xlogx_fast(_wr[r]-1) + xlogx_fast(_wr[nr]+1)) * _mrs.size();
#pragma omp parallel for schedule(runtime) reduction(+:Sa, Sb) \
if (_mrs.size() > OPENMP_MIN_THRESH)
for (size_t i = 0; i < _mrs.size(); ++i)
{
auto& mrsi = _mrs[i];
......@@ -185,9 +191,6 @@ public:
Sb += -2 * (xlogx_fast(mrs) + xlogx_fast(mnrs));
Sa += -2 * (xlogx_fast(mrs-1) + xlogx_fast(mnrs+1));
Sb += xlogx_fast(_wr[r]) + xlogx_fast(_wr[nr]);
Sa += xlogx_fast(_wr[r]-1) + xlogx_fast(_wr[nr]+1);
}
return (Sa - Sb);
}
......@@ -206,7 +209,8 @@ public:
}
// Computes the move proposal probability
double get_move_prob(size_t, size_t r, size_t s, double, double d, bool reverse)
double get_move_prob(size_t, size_t r, size_t s, double, double d,
bool reverse)
{
size_t B = _candidate_blocks.size();
if (reverse)
......@@ -252,11 +256,13 @@ public:
double S = 0, S_n = 0;
for (auto nr : _wr)
S_n += xlogx_fast(nr);
#pragma omp parallel for schedule(runtime) reduction(+:S)
for (size_t i = 0; i < _mrs.size(); ++i)
{
for (auto c : _mrs[i])
for (auto& c : _mrs[i])
S -= 2 * xlogx_fast(c.second);
for (auto rn : _nr[i])
for (auto& rn : _nr[i])
S += xlogx_fast(rn.second);
S += S_n;
}
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2020 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_partition_centroid_rmi.hh"
#include "../support/graph_state.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, RMICenterState, BLOCK_STATE_params)
python::object make_rmi_center_state(boost::python::object ostate)
{
python::object state;
block_state::make_dispatch(ostate,
[&](auto& s){state = python::object(s);});
return state;
}
void export_rmi_center_state()
{
using namespace boost::python;
def("make_rmi_center_state", &make_rmi_center_state);
block_state::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) =
&state_t::virtual_move;
class_<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("log_omega",
+[](python::object onr, python::object ons) -> double
{
auto nr = get_array<int32_t, 1>(onr);
auto ns = get_array<int32_t, 1>(ons);
return log_omega(nr, ns, [](auto& x) {return x;});
});
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2020 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_PARTITION_CENTROID_HH
#define GRAPH_PARTITION_CENTROID_HH
#include "config.h"
#include <vector>
#include "../blockmodel/graph_blockmodel_util.hh"
#include "../support/graph_state.hh"
#include "openmp_lock.hh"
namespace graph_tool
{
using namespace boost;
using namespace std;
template <class Vr, class Vs, class T>
double log_omega(Vr&& nr, Vs&& ns, T&& get_n)
{
double S = 0;
size_t N = 0;
size_t Bx = 0;
size_t By = 0;
for (auto& n : nr)
{
Bx += (get_n(n) > 0);
N += get_n(n);
}
for (auto& n : ns)
By += (get_n(n) > 0);
S += (Bx - 1) * (By - 1) * log(N + (Bx * By) / 2.);
double w = N / (N + (Bx * By) / 2.);
double y2 = 0, lny = 0;
double x2 = 0, lnx = 0;
for (auto& n : nr)
{
double x_r = (1-w)/Bx + (w * get_n(n))/N;
x2 += x_r * x_r;
lnx += log(x_r);
}
for (auto& n : ns)
{
double y_s = (1-w)/By + (w * get_n(n))/N;
y2 += y_s * y_s;
lny += log(y_s);
}
double mu = (Bx + 1) / (Bx * y2) - 1./Bx;
double nu = (By + 1) / (By * x2) - 1./By;
S += (Bx + nu - 2) * lny / 2;
S += (By + mu - 2) * lnx / 2;
S += lgamma(mu * Bx)/2 + lgamma(nu * By)/2;
S -= By * (lgamma(nu) + lgamma(Bx))/2. + Bx * (lgamma(mu) + lgamma(By))/2.;
return S;
}
typedef multi_array_ref<int32_t,2> bs_t;
typedef multi_array_ref<int32_t,1> b_t;
#define BLOCK_STATE_params \
((g, &, always_directed_never_reversed, 1)) \
((_abg, &, boost::any&, 0)) \
((bs,, bs_t, 0)) \
((b,, b_t, 0))
GEN_STATE_BASE(RMICenterStateBase, BLOCK_STATE_params)
template <class... Ts>
class RMICenterState
: public RMICenterStateBase<Ts...>
{
public:
GET_PARAMS_USING(RMICenterStateBase<Ts...>, BLOCK_STATE_params)
GET_PARAMS_TYPEDEF(Ts, BLOCK_STATE_params)
template <class... ATs,
typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
RMICenterState(ATs&&... args)
: RMICenterStateBase<Ts...>(std::forward<ATs>(args)...),
_bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
_mrs(_bs.shape()[0]),
_nr(_bs.shape()[0]),
_N(_bs.shape()[1]),
_wr(_N),
_empty_pos(_N),
_candidate_pos(_N),
_bclabel(_N),
_pclabel(_N)
{
for (size_t r : _b)
_wr[r]++;
for (size_t r = 0; r < _N; ++r)
{
if (_wr[r] == 0)
add_element(_empty_blocks, _empty_pos, r);
else
add_element(_candidate_blocks, _candidate_pos, r);
}
for (size_t i = 0; i < _mrs.size(); ++i)
{
for (size_t v = 0; v < _N; ++v)
{
auto r = _b[v];
auto s = _bs[i][v];
_mrs[i][{r,s}]++;
_nr[i][s]++;
}
}
}
typedef typename
std::conditional<is_directed_::apply<g_t>::type::value,
GraphInterface::multigraph_t,
undirected_adaptor<GraphInterface::multigraph_t>>::type
bg_t;
bg_t& _bg;
std::vector<gt_hash_map<std::tuple<size_t, size_t>, size_t>> _mrs;
std::vector<gt_hash_map<size_t, size_t>> _nr;
size_t _N;
std::vector<size_t> _wr;
std::vector<size_t> _empty_blocks;
std::vector<size_t> _empty_pos;
std::vector<size_t> _candidate_blocks;
std::vector<size_t> _candidate_pos;
std::vector<size_t> _bclabel;
std::vector<size_t> _pclabel;
constexpr static BlockStateVirtualBase* _coupled_state = nullptr;
typedef int m_entries_t;
bool _egroups_update = true;
// =========================================================================
// State modification
// =========================================================================
void move_vertex(size_t v, size_t nr)
{
size_t r = _b[v];
if (nr == r)
return;
_wr[r]--;
_wr[nr]++;
#pragma omp parallel for schedule(runtime) \
if (_mrs.size() > OPENMP_MIN_THRESH)
for (size_t i = 0; i < _mrs.size(); ++i)
{
auto& mrsi = _mrs[i];
size_t s = _bs[i][v];
auto iter = mrsi.find({r, s});
assert(iter != mrsi.end());
iter->second--;
if (iter->second == 0)
mrsi.erase(iter);
mrsi[{nr,s}]++;
}
if (_wr[r] == 0)
{
add_element(_empty_blocks, _empty_pos, r);
remove_element(_candidate_blocks, _candidate_pos, r);
}
if (_wr[nr] == 1)
{
remove_element(_empty_blocks, _empty_pos, nr);
add_element(_candidate_blocks, _candidate_pos, nr);
}
_b[v] = nr;
}
size_t virtual_remove_size(size_t v)
{
return _wr[_b[v]] - 1;
}
constexpr void add_block(size_t)
{
}
double virtual_move(size_t v, size_t r, size_t nr)
{
if (r == nr)
return 0;
double Sb = entropy();
move_vertex(v, nr);
double Sa = entropy();
move_vertex(v, r);
return (Sa - Sb);
}
size_t get_empty_block(size_t, bool)
{
return _empty_blocks.back();
}
size_t sample_block(size_t, double, double d, rng_t& rng)
{
std::bernoulli_distribution new_r(d);
if (d > 0 && !_empty_blocks.empty() && new_r(rng))
return uniform_sample(_empty_blocks, rng);
return uniform_sample(_candidate_blocks, rng);
}
// Computes the move proposal probability
double get_move_prob(size_t, size_t r, size_t s, double, double d,
bool reverse)
{
size_t B = _candidate_blocks.size();
if (reverse)
{
if (_wr[s] == 1)
return d;
if (_wr[r] == 0)
B++;
}
else
{
if (_wr[s] == 0)
return d;
}
if (B == _N)
d = 0;
return (1. - d) / B;
}
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse, MEntries&&)
{
return get_move_prob(v, r, s, c, d, reverse);
}
template <class EArgs>
double virtual_move(size_t v, size_t r, size_t nr, EArgs&&)
{
return virtual_move(v, r, nr);
}
template <class EArgs, class MEntries>
double virtual_move(size_t v, size_t r, size_t nr, EArgs&&, MEntries&&)
{
return virtual_move(v, r, nr);
}
double entropy()
{
double S = 0, S_n = 0;
gt_hash_map<int, int> wr;
for (auto r : _candidate_blocks)
{
if (_wr[r] == 0)
continue;
S_n -= lgamma_fast(_wr[r] + 1);
wr[r] = _wr[r];
}
S += (lgamma_fast(_N + 1) + S_n) * _mrs.size();
#pragma omp parallel for schedule(runtime) reduction(+:S)
for (size_t i = 0; i < _mrs.size(); ++i)
{
for (auto& c : _mrs[i])
S += lgamma_fast(c.second + 1);
for (auto& rn : _nr[i])
S -= lgamma_fast(rn.second + 1);
S -= log_omega(wr, _nr[i], [](auto& x) {return x.second;});
}
return -S;
}
void init_mcmc(double, double)
{
}
constexpr size_t node_weight(size_t)
{
return 1;
}
bool is_last(size_t v)
{
return _wr[_b[v]] == 1;
}
constexpr bool allow_move(size_t, size_t)
{
return true;
}
template <class V>
void push_state(V&) {}
void pop_state() {}
void store_next_state(size_t) {}
void clear_next_state() {}
};
} // graph_tool namespace
#endif //GRAPH_PARTITION_CENTROID_HH
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2020 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_partition_centroid_rmi.hh"
#include "../blockmodel/graph_blockmodel_mcmc.hh"
#include "../loops/mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, RMICenterState, BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object rmi_mcmc_sweep(python::object omcmc_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;