Commit 694f4369 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

uncertain SBM: fix edge sampling

parent 79dd3b1f
......@@ -6,6 +6,7 @@ from pylab import *
from graph_tool.all import *
import numpy.random
from numpy.random import randint
from collections import defaultdict
import scipy.stats
import os.path
......@@ -115,6 +116,46 @@ for directed in [True, False]:
savefig("test_mcmc/test_mcmc_move_prob_directed%s.pdf" % directed)
for directed in [True, False]:
for deg_corr in [True, False]:
print(f"test edge sampling directed={directed}, deg_corr={deg_corr}...")
g = collection.data["football"]
g.set_directed(directed)
state = LatentMultigraphBlockState(g, state_args=dict(deg_corr=deg_corr))
for i in range(1000):
state.multiflip_mcmc_sweep(niter=10)
print(state)
esampler = state.bstate._state.get_edge_sampler(False)
ecount = defaultdict(int)
for i in range(10000000):
u, v = esampler.sample(graph_tool._get_rng())
if v < u and not directed:
u, v = v, u
ecount[(u,v)] += 1
N = sum([x for x in ecount.values()])
Lp = []
Ls = []
for e, c in ecount.items():
p = log(c)-log(N)
u, v = e
e = state.bstate.g.edge(u, v)
if e is None:
m = 0
else:
m = state.bstate.eweight[e]
Le = esampler.log_prob(u, v, m)
Lp.append(Le)
Ls.append(p)
clf()
plot(Ls, Lp, "o")
plot(Ls, Ls, "-")
savefig(f"test_mcmc/test_mcmc_edge_sample_directed{directed}_deg_corr{deg_corr}.pdf")
if os.path.exists("g_small.gt"):
g_small = load_graph("g_small.gt")
else:
......
......@@ -22,6 +22,7 @@
#include "graph_blockmodel_util.hh"
#include "graph_blockmodel.hh"
#include "../uncertain/graph_blockmodel_sample_edge.hh"
#include "../support/graph_state.hh"
using namespace boost;
......@@ -98,6 +99,27 @@ void export_sbm_state()
.def("rebuild_neighbor_sampler",
&state_t::rebuild_neighbor_sampler)
.def("sync_emat",
&state_t::sync_emat);
&state_t::sync_emat)
.def("get_edge_sampler",
+[](state_t& state, bool edges_only)
{
return SBMEdgeSampler(state, edges_only);
});
typedef SBMEdgeSampler<state_t> edge_sampler_t;
class_<edge_sampler_t>
es(name_demangle(typeid(edge_sampler_t).name()).c_str(),
no_init);
es.def("sample",
+[](edge_sampler_t& es, rng_t& rng)
{
auto e = es.sample(rng);
return python::make_tuple(get<0>(e), get<1>(e));
})
.def("log_prob",
+[](edge_sampler_t& es, size_t u, size_t v, size_t m)
{
return es.log_prob(u, v, m, 0);
});
});
}
......@@ -236,7 +236,7 @@ struct MCMC
double a = 0;
if (dm != 0)
a += (_edge_sampler.log_prob(u, v, m + dm, dm) -
a += (_edge_sampler.log_prob(u, v, m, dm) -
_edge_sampler.log_prob(u, v, m, 0));
a -= log(get_move_prob(dm, dx, m, x));
......@@ -256,17 +256,14 @@ struct MCMC
{
_state.update_edge(u, v, dx);
}
else if (dm < 0)
{
size_t m = get<0>(node_state(u, v));
_edge_sampler.template update_edge<false>(u, v, m);
_state.remove_edge(u, v);
}
else
{
size_t m = get<0>(node_state(u, v));
_state.add_edge(u, v, dx);
_edge_sampler.template update_edge<true>(u, v, m);
_edge_sampler.update_edge(u, v, m, dm);
if (dm < 0)
_state.remove_edge(u, v);
else
_state.add_edge(u, v, dx);
}
}
......
......@@ -102,6 +102,9 @@ public:
_NB = B * B;
}
SBMEdgeSampler(const SBMEdgeSampler& other)
: SBMEdgeSampler(other._state, other._edges_only) {}
std::tuple<size_t, size_t> get_edge(size_t u, size_t v)
{
if (!graph_tool::is_directed(_state._g) && u > v)
......@@ -110,38 +113,31 @@ public:
}
template <bool add>
void update_edge(size_t u, size_t v, size_t m)
void update_edge(size_t u, size_t v, size_t m, int delta)
{
if (_edges_only)
return;
if (add)
if (m == 0 && delta > 0)
{
if (m == 0)
{
_edges.push_back(get_edge(u, v));
_edge_pos[_edges.back()] = _edges.size() - 1;
}
_edges.push_back(get_edge(u, v));
_edge_pos[_edges.back()] = _edges.size() - 1;
}
else
if (m > 0 && m + delta == 0)
{
if (m == 1)
auto iter = _edge_pos.find(get_edge(u, v));
size_t pos = iter->second;
_edge_pos.erase(iter);
if (pos < _edges.size() - 1)
{
auto iter = _edge_pos.find(get_edge(u, v));
size_t pos = iter->second;
_edge_pos.erase(iter);
if (pos < _edges.size() - 1)
{
std::swap(_edges[pos], _edges.back());
_edge_pos[_edges[pos]] = pos;
}
_edges.pop_back();
std::swap(_edges[pos], _edges.back());
_edge_pos[_edges[pos]] = pos;
}
_edges.pop_back();
}
constexpr int delta = (add) ? 0 : -1;
_E += (add) ? 1 : -1;
_E += delta;
size_t r = _state._b[u];
size_t s = _state._b[v];
......@@ -149,12 +145,15 @@ public:
if (me != _state._emat.get_null_edge())
{
auto ers = _state._mrs[me] + delta;
if (!add || ers > 1)
if (ers == 0)
{
_rs_sampler.remove(_rs_pos[me]);
if (ers > 0)
_rs_pos[me] = _rs_sampler.insert({r,s}, ers);
else
_rs_pos[me] = std::numeric_limits<size_t>::max();
}
else
{
_rs_pos[me] = _rs_sampler.insert({r,s}, ers);
}
}
if (_state._deg_corr)
......@@ -207,9 +206,6 @@ public:
template <class RNG>
std::tuple<size_t, size_t> sample(RNG& rng)
{
// std::uniform_int_distribution<size_t> sample(0, _N-1);
// return {sample(rng), sample(rng)};
if (_edges_only)
{
std::bernoulli_distribution coin(_edges.size() /
......@@ -329,11 +325,11 @@ public:
if (!graph_tool::is_directed(g) && u != v)
lp += log(2);
if (m > 0)
if (m + delta > 0)
{
double rp;
if (m == 1 && delta == 1)
rp = -log(_edges.size() + delta);
if (m == 0)
rp = -log(_edges.size() + 1);
else
rp = -log(_edges.size());
double x = std::max(rp, lp);
......@@ -342,7 +338,7 @@ public:
}
else
{
return lp;
return lp - log(2);
}
}
......
......@@ -105,42 +105,54 @@ struct MCMC
{
_e = _edge_sampler.sample(rng);
std::bernoulli_distribution coin(.5);
size_t m = node_state(get<0>(_e), get<1>(_e));
if (m > 0 && coin(rng))
{
return -1;
}
else
{
return 1;
}
int m = node_state(get<0>(_e), get<1>(_e));
std::geometric_distribution<int> sample_m(1./(m + 2));
return sample_m(rng) - m;
}
std::tuple<double, double>
virtual_move_dS(size_t, int dm)
{
if (dm == 0)
return {0., 0.};
size_t u, v;
std::tie(u, v) = get_edge();
double dS = 0;
if (dm < 0)
{
dS = _state.remove_edge_dS(u, v, _entropy_args);
for (int i = 0; i < -dm-1; ++i)
{
_state.remove_edge(u, v);
dS += _state.remove_edge_dS(u, v, _entropy_args);
}
for (int i = 0; i < -dm-1; ++i)
_state.add_edge(u, v);
}
else
{
dS = _state.add_edge_dS(u, v, _entropy_args);
for (int i = 0; i < dm-1; ++i)
{
_state.add_edge(u, v);
dS += _state.add_edge_dS(u, v, _entropy_args);
}
for (int i = 0; i < dm-1; ++i)
_state.remove_edge(u, v);
}
size_t m = node_state(u, v);
double a = (_edge_sampler.log_prob(u, v, m + dm, dm) -
double a = (_edge_sampler.log_prob(u, v, m, dm) -
_edge_sampler.log_prob(u, v, m, 0));
if (m > 0)
{
a -= -log(2);
}
if (m + dm > 0)
{
a += -log(2);
}
size_t nm = m + dm;
a -= nm * safelog_fast(m + 1) - (nm + 1) * safelog_fast(m + 2);
a += m * safelog_fast(nm + 1) - (m + 1) * safelog_fast(nm + 2);
return std::make_tuple(dS, a);
}
......@@ -152,15 +164,18 @@ struct MCMC
size_t u, v;
std::tie(u, v) = get_edge();
size_t m = node_state(u, v);
_edge_sampler.update_edge(u, v, m, dm);
if (dm < 0)
{
_edge_sampler.template update_edge<false>(u, v, m);
_state.remove_edge(u, v);
for (int i = 0; i < -dm; ++i)
_state.remove_edge(u, v);
}
else
{
_state.add_edge(u, v);
_edge_sampler.template update_edge<true>(u, v, m);
for (int i = 0; i < dm; ++i)
_state.add_edge(u, v);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment