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 0cebc67b authored by Tiago Peixoto's avatar Tiago Peixoto

inference.blockmodel: Implement multiple edge covariates

parent 5053de1f
......@@ -475,6 +475,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_core)
boost::mpl::for_each<boost::mpl::push_back<scalar_types,string>::type>(export_vector_types());
export_vector_types()(size_t(), "size_t");
export_vector_types()(std::vector<double>(), "Vector_double");
class_<GraphInterface>("GraphInterface", init<>())
.def(init<GraphInterface,bool,boost::python::object,
......
......@@ -120,7 +120,8 @@ struct in_degreeS
auto get_in_degree(typename boost::graph_traits<Graph>::vertex_descriptor v,
const Graph& g, std::true_type, Weight& weight) const
{
typename boost::property_traits<Weight>::value_type d = 0;
typedef typename boost::property_traits<Weight>::value_type val_t;
val_t d = val_t();
typename boost::graph_traits<Graph>::in_edge_iterator e, e_end;
for (std::tie(e, e_end) = in_edges(v, g); e != e_end; ++e)
d += get(weight, *e);
......@@ -179,7 +180,8 @@ struct out_degreeS
auto get_out_degree(typename boost::graph_traits<Graph>::vertex_descriptor v,
const Graph& g, const Weight& weight) const
{
typename boost::property_traits<Weight>::value_type d = 0;
typedef typename boost::property_traits<Weight>::value_type val_t;
val_t d = val_t();
typename boost::graph_traits<Graph>::out_edge_iterator e, e_end;
for (std::tie(e, e_end) = out_edges(v, g); e != e_end; ++e)
d += get(weight, *e);
......
......@@ -67,18 +67,17 @@ void clear_xlogx()
void init_lgamma(size_t x)
{
using namespace boost::math::policies;
#pragma omp critical (_lgamma_)
{
size_t old_size = __lgamma_cache.size();
if (x >= old_size)
{
__lgamma_cache.resize(x + 1);
if (old_size == 0)
__lgamma_cache[0] = numeric_limits<double>::infinity();
for (size_t i = std::max(old_size, size_t(1)); i < __lgamma_cache.size(); ++i)
__lgamma_cache[i] = boost::math::lgamma(i);
for (size_t i = std::max(old_size, size_t(1));
i < __lgamma_cache.size(); ++i)
__lgamma_cache[i] = lgamma(i);
}
}
}
......
......@@ -220,6 +220,7 @@ void export_blockmodel_state()
.def_readwrite("dense", &entropy_args_t::dense)
.def_readwrite("multigraph", &entropy_args_t::multigraph)
.def_readwrite("adjacency", &entropy_args_t::adjacency)
.def_readwrite("recs", &entropy_args_t::recs)
.def_readwrite("partition_dl", &entropy_args_t::partition_dl)
.def_readwrite("degree_dl", &entropy_args_t::degree_dl)
.def_readwrite("degree_dl_kind", &entropy_args_t::degree_dl_kind)
......@@ -232,10 +233,11 @@ void export_blockmodel_state()
enum_<weight_type>("rec_type")
.value("none", weight_type::NONE)
.value("positive", weight_type::POSITIVE)
.value("signed", weight_type::SIGNED)
.value("real_exponential", weight_type::REAL_EXPONENTIAL)
.value("real_normal", weight_type::REAL_NORMAL)
.value("discrete_geometric", weight_type::DISCRETE_GEOMETRIC)
.value("discrete_poisson", weight_type::DISCRETE_POISSON)
.value("discrete_binomial", weight_type::DISCRETE_BINOMIAL)
.value("delta_t", weight_type::DELTA_T);
def("make_block_state", &make_block_state);
......@@ -255,4 +257,10 @@ void export_blockmodel_state()
def("log_q_approx", log_q_approx);
def("log_q_approx_big", log_q_approx_big);
def("log_q_approx_small", log_q_approx_small);
def("positive_w_log_P", positive_w_log_P<size_t>);
def("signed_w_log_P", signed_w_log_P<size_t>);
def("geometric_w_log_P", geometric_w_log_P<size_t>);
def("binomial_w_log_P", binomial_w_log_P<size_t>);
def("poisson_w_log_P", poisson_w_log_P<size_t>);
}
This diff is collapsed.
......@@ -462,12 +462,13 @@ struct Layers
}
double entropy(bool dense, bool multigraph, bool deg_entropy,
bool exact)
bool exact, bool recs)
{
double S = 0;
if (_master)
{
S += BaseState::entropy(dense, multigraph, deg_entropy, exact);
S += BaseState::entropy(dense, multigraph, deg_entropy, exact,
recs);
S -= covariate_entropy(_bg, _mrs);
if (multigraph)
S -= BaseState::get_parallel_entropy();
......@@ -481,7 +482,8 @@ struct Layers
else
{
for (auto& state : _layers)
S += state.entropy(dense, multigraph, deg_entropy, exact);
S += state.entropy(dense, multigraph, deg_entropy, exact,
recs);
}
return S;
}
......
......@@ -1248,12 +1248,17 @@ public:
void set_move(size_t, size_t) {}
template <class... DVals>
void insert_delta(size_t t, size_t s, DVals... delta)
void insert_delta(bool add, size_t t, size_t s, DVals... delta)
{
if (!is_directed::apply<Graph>::type::value && (t > s))
std::swap(t, s);
_entries[_pos] = make_pair(t, s);
add_to_tuple(_delta[_pos], delta...);
if (add)
tuple_op(_delta[_pos], [&](auto& r, auto& v){ r += v; },
delta...);
else
tuple_op(_delta[_pos], [&](auto& r, auto& v){ r -= v; },
delta...);
++_pos;
}
......@@ -1314,6 +1319,8 @@ public:
return emat.get_me(t, s);
}
std::tuple<EVals...> _self_weight;
private:
size_t _pos;
std::array<pair<size_t, size_t>, 2> _entries;
......
......@@ -39,26 +39,61 @@
#include <omp.h>
#endif
namespace std
{
template <class T>
void operator+=(vector<T>& ret, const vector<T>& v)
{
ret.resize(max(ret.size(), v.size()));
for (size_t i = 0; i < v.size(); ++i)
ret[i] += v[i];
}
template <class T>
void operator-=(vector<T>& ret, const vector<T>& v)
{
ret.resize(max(ret.size(), v.size()));
for (size_t i = 0; i < v.size(); ++i)
ret[i] -= v[i];
}
template <class T1, class T2>
void operator*=(vector<T1>& ret, const T2& v)
{
for (auto& x : ret)
x *= v;
}
template <class T1, class T2>
void operator/=(vector<T1>& ret, const T2& v)
{
for (auto& x : ret)
x /= v;
}
}
namespace graph_tool
{
// tuple utils
template <size_t i, class T>
void add_to_tuple_imp(T&)
template <size_t i, class T, class OP>
void tuple_op_imp(T&, OP&&)
{
}
template <size_t i, class T, class Ti, class... Ts>
void add_to_tuple_imp(T& tuple, Ti&& v, Ts&&... vals)
template <size_t i, class T, class OP, class Ti, class... Ts>
void tuple_op_imp(T& tuple, OP&& op, Ti&& v, Ts&&... vals)
{
get<i>(tuple) += v;
add_to_tuple_imp<i+1>(tuple, vals...);
op(get<i>(tuple), v);
tuple_op_imp<i+1>(tuple, op, vals...);
}
template <class T, class... Ts>
void add_to_tuple(T& tuple, Ts&&... vals)
template <class OP, class T, class... Ts>
__attribute__((flatten))
void tuple_op(T& tuple, OP&& op, Ts&&... vals)
{
add_to_tuple_imp<0>(tuple, vals...);
tuple_op_imp<0>(tuple, op, vals...);
}
namespace detail {
......@@ -95,6 +130,7 @@ struct entropy_args_t
bool multigraph;
bool exact;
bool adjacency;
bool recs;
bool partition_dl;
bool degree_dl;
deg_dl_kind degree_dl_kind;
......@@ -222,8 +258,8 @@ double positive_w_log_P(DT N, double x, double alpha, double beta)
{
if (N == 0)
return 0.;
return boost::math::lgamma(N + alpha) - boost::math::lgamma(alpha)
+ alpha * log(beta) - (alpha + N) * log(beta + x);
return lgamma(N + alpha) - lgamma(alpha) + alpha * log(beta) -
(alpha + N) * log(beta + x);
}
// normal
......@@ -236,9 +272,9 @@ double signed_w_log_P(DT N, double x, double v, double m0, double k0, double v0,
auto k_n = k0 + N;
auto nu_n = nu0 + N;
auto v_n = (v0 * nu0 + v + ((N * k0)/(k0 + N)) * pow(m0 - x/N, 2.)) / nu_n;
return boost::math::lgamma(nu_n/2.) - boost::math::lgamma(nu0/2.)
+ (log(k0) - log(k_n))/2. + (nu0 / 2.) * log(nu0 * v0)
- (nu_n / 2.) * log(nu_n * v_n) - (N/2.) * log(M_PI);
return lgamma(nu_n / 2.) - lgamma(nu0 / 2.) + (log(k0) - log(k_n)) / 2. +
(nu0 / 2.) * log(nu0 * v0) - (nu_n / 2.) * log(nu_n * v_n) - (N / 2.) *
log(M_PI);
}
// discrete: geometric
......@@ -247,21 +283,26 @@ double geometric_w_log_P(DT N, double x, double alpha, double beta)
{
if (N == 0)
return 0.;
return boost::math::lgamma(alpha + beta)
+ boost::math::lgamma(alpha + 1) + boost::math::lgamma(beta + 1)
- boost::math::lgamma(alpha) - boost::math::lgamma(beta)
- boost::math::lgamma(x + alpha + beta + 1);
return lbeta(N + alpha, x + beta) - lbeta(alpha, beta);
}
// discrete: binomial
template <class DT>
double binomial_w_log_P(DT N, double x, size_t n, double alpha, double beta)
{
if (N == 0)
return 0.;
return lbeta(x + alpha, N * n - x + beta) - lbeta(alpha, beta);
}
// discrete: Poisson
template <class DT>
double poisson_w_log_P(DT N, double x, double alpha, double beta)
double poisson_w_log_P(DT N, double x, double r, double theta)
{
if (N == 0)
return 0.;
return boost::math::lgamma(x + alpha)
- boost::math::lgamma(x + 1) - boost::math::lgamma(alpha)
+ alpha * log(beta) - (x + alpha) * log(beta + 1);
return lgamma(x + r) - lgamma(r) - r * log(theta) - (x + r) *
log(N + 1. / theta);
}
// ===============
......@@ -1012,14 +1053,15 @@ public:
template <class... DVals>
__attribute__((flatten))
void insert_delta(size_t t, size_t s, DVals... delta)
void insert_delta(bool add, size_t t, size_t s, DVals&&... delta)
{
insert_delta_imp(t, s, typename is_directed::apply<Graph>::type(),
delta...);
insert_delta_imp(add, t, s, typename is_directed::apply<Graph>::type(),
std::forward<DVals>(delta)...);
}
template <class... DVals>
void insert_delta_imp(size_t t, size_t s, std::true_type, DVals... delta)
void insert_delta_imp(bool add, size_t t, size_t s, std::true_type,
DVals&&... delta)
{
bool src = false;
if (t != _rnr.first && t != _rnr.second)
......@@ -1043,11 +1085,17 @@ public:
_entries.emplace_back(t, s);
_delta.emplace_back();
}
add_to_tuple(_delta[field[s]], delta...);
if (add)
tuple_op(_delta[field[s]], [&](auto& r, auto& v){ r += v; },
delta...);
else
tuple_op(_delta[field[s]], [&](auto& r, auto& v){ r -= v; },
delta...);
}
template <class... DVals>
void insert_delta_imp(size_t t, size_t s, std::false_type, DVals... delta)
void insert_delta_imp(bool add, size_t t, size_t s, std::false_type,
DVals&&... delta)
{
if (t > s)
std::swap(t, s);
......@@ -1066,7 +1114,12 @@ public:
_entries.emplace_back(t, s);
_delta.emplace_back();
}
add_to_tuple(_delta[field[s]], delta...);
if (add)
tuple_op(_delta[field[s]], [&](auto& r, auto& v){ r += v; },
delta...);
else
tuple_op(_delta[field[s]], [&](auto& r, auto& v){ r -= v; },
delta...);
}
size_t get_field(size_t r, size_t s)
......@@ -1152,6 +1205,8 @@ public:
return _mes[field];
}
std::tuple<EVals...> _self_weight;
private:
static constexpr size_t _null = numeric_limits<size_t>::max();
static const std::tuple<EVals...> _null_delta;
......@@ -1184,8 +1239,17 @@ void modify_entries(Vertex v, Vertex r, Vprop& _b, Graph& g, Eprop& eweights,
Eprops&... eprops)
{
typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
std::tuple<int, typename property_traits<Eprops>::value_type...>
self_weight;
auto& self_weight = m_entries._self_weight;
if (!is_directed::apply<Graph>::type::value)
{
tuple_apply([&](auto&&... vals)
{
auto op = [](auto& x) -> auto& { x *= 0; return x; };
auto f = [](auto&...) {};
f(op(vals)...);
}, self_weight);
}
for (auto e : out_edges_range(v, g))
{
if (efilt(e))
......@@ -1198,26 +1262,22 @@ void modify_entries(Vertex v, Vertex r, Vprop& _b, Graph& g, Eprop& eweights,
if (Add && u == v)
s = r;
if (Add)
m_entries.insert_delta(r, s, ew, eprops[e]...);
else
m_entries.insert_delta(r, s, -ew, -eprops[e]...);
m_entries.insert_delta(Add, r, s, ew, eprops[e]...);
if ((u == v || is_loop(v)) && !is_directed::apply<Graph>::type::value)
add_to_tuple(self_weight, ew, eprops[e]...);
tuple_op(self_weight, [&](auto& r, auto& v){ r += v; },
ew, eprops[e]...);
}
if (get<0>(self_weight) > 0 && get<0>(self_weight) % 2 == 0 &&
!is_directed::apply<Graph>::type::value)
{
tuple_apply([&](auto&&... vals)
{
if (Add)
m_entries.insert_delta(r, r,
(-vals / 2)...);
else
m_entries.insert_delta(r, r,
(vals / 2)...);
auto op = [](auto& x) -> auto& { x /= 2; return x; };
m_entries.insert_delta(!Add, r, r, op(vals)...);
}, self_weight);
}
for (auto e : in_edges_range(v, g))
{
......@@ -1229,10 +1289,7 @@ void modify_entries(Vertex v, Vertex r, Vprop& _b, Graph& g, Eprop& eweights,
vertex_t s = _b[u];
int ew = eweights[e];
if (Add)
m_entries.insert_delta(s, r, ew, eprops[e]...);
else
m_entries.insert_delta(s, r, -ew, -eprops[e]...);
m_entries.insert_delta(Add, s, r, ew, eprops[e]...);
}
}
......
......@@ -28,7 +28,6 @@ using namespace graph_tool;
template <class Value>
void vector_map(boost::python::object ovals, boost::python::object omap)
{
multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
multi_array_ref<Value,1> map = get_array<Value,1>(omap);
......@@ -45,7 +44,6 @@ void vector_map(boost::python::object ovals, boost::python::object omap)
template <class Value>
void vector_continuous_map(boost::python::object ovals)
{
multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
gt_hash_map<Value, size_t> map;
......@@ -62,7 +60,6 @@ void vector_continuous_map(boost::python::object ovals)
template <class Value>
void vector_rmap(boost::python::object ovals, boost::python::object omap)
{
multi_array_ref<Value,1> vals = get_array<Value,1>(ovals);
multi_array_ref<Value,1> map = get_array<Value,1>(omap);
......
......@@ -29,33 +29,28 @@ namespace graph_tool
{
using namespace boost;
// Warning: std::lgamma(x) is not thread-safe! Although in the context of this
// program the outputs should _always_ be positive, we use boost::math::lgamma
// instead.
inline double lbinom(double N, double k)
{
if (N == 0 || k == 0 || k > N)
return 0;
assert(N > 0);
assert(k > 0);
return ((boost::math::lgamma(N + 1) - boost::math::lgamma(k + 1))
- boost::math::lgamma(N - k + 1));
return ((lgamma(N + 1) - lgamma(k + 1)) - lgamma(N - k + 1));
}
inline double lbinom_fast(int N, int k)
{
if (N == 0 || k == 0 || k > N)
return 0;
return lgamma_fast(N + 1) - lgamma_fast(N - k + 1) - lgamma_fast(k + 1);
return ((lgamma_fast(N + 1) - lgamma_fast(k + 1)) - lgamma_fast(N - k + 1));
}
inline double lbinom_careful(double N, double k)
{
if (N == 0 || k == 0 || k >= N)
return 0;
double lgN = boost::math::lgamma(N + 1);
double lgk = boost::math::lgamma(k + 1);
double lgN = lgamma(N + 1);
double lgk = lgamma(k + 1);
if (lgN - lgk > 1e8)
{
// We have N >> k. Use Stirling's approximation: ln N! ~ N ln N - N
......@@ -64,10 +59,15 @@ inline double lbinom_careful(double N, double k)
}
else
{
return lgN - boost::math::lgamma(N - k + 1) - lgk;
return lgN - lgamma(N - k + 1) - lgk;
}
}
inline double lbeta(double x, double y)
{
return lgamma(x) + lgamma(y) - lgamma(x + y);
}
template <class Vec, class PosMap, class Val>
void remove_element(Vec& vec, PosMap& pos, Val val)
{
......
This diff is collapsed.
......@@ -31,6 +31,7 @@ import numpy
from collections import defaultdict
from scipy.special import gammaln
import copy
import warnings
from .. import group_vector_property, ungroup_vector_property, Vector_size_t, \
perfect_prop_hash
......@@ -283,8 +284,8 @@ class LayeredBlockState(OverlapBlockState, BlockState):
kwargs.pop("rec_params", None)
if len(kwargs) > 0:
raise ValueError("unrecognized keyword arguments: " +
str(list(kwargs.keys())))
warnings.warn("unrecognized keyword arguments: " +
str(list(kwargs.keys())))
def get_N(self):
r"Returns the total number of edges."
......
......@@ -98,7 +98,8 @@ class NestedBlockState(object):
degree_dl=True,
degree_dl_kind="distributed",
edges_dl=True,
exact=True),
exact=True,
recs=False),
**hentropy_args)
self.levels = [base_type(g, b=bs[0], **self.state_args)]
for i, b in enumerate(bs[1:]):
......
......@@ -271,8 +271,8 @@ class OverlapBlockState(BlockState):
exact=True)
if len(kwargs) > 0:
raise ValueError("unrecognized keyword arguments: " +
str(list(kwargs.keys())))
warnings.warn("unrecognized keyword arguments: " +
str(list(kwargs.keys())))
def __repr__(self):
return "<OverlapBlockState object with %d blocks,%s for graph %s, at 0x%x>" % \
......@@ -351,15 +351,23 @@ class OverlapBlockState(BlockState):
vweight=wr if vweight else None,
b=bg.vertex_index.copy("int") if b is None else b,
deg_corr=deg_corr,
rec_type=kwargs.get("rec_type", self.rec_type if vweight else None),
rec=kwargs.get("rec", bg.ep.rec if vweight else None),
drec=kwargs.get("drec", bg.ep.drec if vweight else None),
rec_params=kwargs.get("rec_params", self.rec_params),
allow_empty=kwargs.get("allow_empty",
rec_types=kwargs.pop("rec_types",
self.rec_types if vweight else None),
recs=kwargs.pop("recs",
ungroup_vector_property(bg.ep.rec,
range(len(self.rec_types)))
if (vweight and
len(self.rec_types) > 0)
else None),
drec=kwargs.pop("drec",
bg.ep.drec if (vweight and
len(self.rec_types) > 0)
else None),
rec_params=kwargs.pop("rec_params", self.rec_params),
allow_empty=kwargs.pop("allow_empty",
self.allow_empty),
max_BE=self.max_BE,
**dmask(kwargs, ["allow_empty", "rec_type",
"rec_params", "rec", "drec"]))
**kwargs)
return state
def get_E(self):
......
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