Commit 0e7a05fb authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Fix numerical issue in OverlapBlockState.entropy() with excessive overlaps

parent bc21f238
......@@ -757,6 +757,9 @@ void export_blockmodel()
def("get_xi_fast", get_xi_fast<double,double>);
def("get_mu_l", python_get_mu_l);
def("polylog", polylog<double>);
def("lbinom_careful", lbinom_careful);
def("lbinom_fast", lbinom_fast);
def("lbinom", lbinom);
def("get_vector", get_vector);
def("vector_map", vector_map<int32_t>);
......
......@@ -284,12 +284,15 @@ struct entropy_parallel_edges
// Dense entropy
// =============
// Warning: lgamma(x) is not thread-safe! However, since in the context of this
// program the outputs should _always_ be positive, this can be overlooked.
__attribute__((always_inline))
inline double lbinom(double N, double k)
{
if (N == 0 || k == 0 || k > N)
return 0;
return lgamma(N + 1) - lgamma(N - k + 1) - lgamma(k + 1);
return (lgamma(N + 1) - lgamma(k + 1)) - lgamma(N - k + 1);
}
__attribute__((always_inline))
......@@ -300,6 +303,25 @@ inline double lbinom_fast(int N, int k)
return lgamma_fast(N + 1) - lgamma_fast(N - k + 1) - lgamma_fast(k + 1);
}
__attribute__((always_inline))
inline double lbinom_careful(double N, double k)
{
if (N == 0 || k == 0 || k >= N)
return 0;
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
// and reorder
return - N * log1p(-k / N) - k * log1p(-k / N) - k - lgk + k * log(N);
}
else
{
return lgN - lgamma(N - k + 1) - lgk;
}
}
// "edge" term of the entropy
template <class Graph>
......
......@@ -560,7 +560,7 @@ struct overlap_partition_stats_t
if (nd == 0)
continue;
double x = lbinom_fast(_actual_B, d);
double ss = lbinom((exp(x) + nd) - 1, nd); // not fast
double ss = lbinom_careful((exp(x) + nd) - 1, nd); // not fast
if (std::isinf(ss) || std::isnan(ss))
ss = nd * x - lgamma_fast(nd + 1);
assert(!std::isinf(ss));
......
......@@ -717,6 +717,12 @@ def model_entropy(B, N, E, directed=False, nr=None):
def lbinom(n, k):
return scipy.special.gammaln(float(n + 1)) - scipy.special.gammaln(float(n - k + 1)) - scipy.special.gammaln(float(k + 1))
def lbinom_careful(n, k):
return libcommunity.lbinom_careful(n, k)
def lbinom_fast(n, k):
return libcommunity.lbinom_fast(n, k)
def partition_entropy(B, N, nr=None):
if nr is None:
S = N * log(B) + log1p(-(1 - 1./B) ** N)
......
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