# Source code for biSBM.int_part

"""
The int_part module computes the number of restricted integer partitions :math:q(m, n) of the integer :math:m
into at most :math:n parts.

Specifically, it counts the number of different degree counts with the sum of degrees
being exactly :math:m and that have at most :math:n nonzero counts. Since the quantity can only computed via a
recurrence, we pre-compute the values to fill up a memoization table when the number of edges and nodes is not too
large; otherwise, we use accurate asymptotic expressions to efficiently compute the values for large arguments.
[peixoto-nonparametric-2017]_

References
----------
.. [peixoto-nonparametric-2017] Tiago P. Peixoto, "Nonparametric
Bayesian inference of the microcanonical stochastic block model",
Phys. Rev. E 95 012317 (2017), :doi:10.1103/PhysRevE.95.012317,
:arxiv:1610.02703

"""
import numpy as np
from numba import njit

from scipy.special import gammaln, spence, loggamma

# for computing the number of restricted partitions of the integer m into at most n pairs
[docs]def log_q(n, k, __q_cache): """log_q Parameters ---------- n : int k : int __q_cache : :class:numpy.ndarray Returns ------- """ n = int(n) k = int(k) if n <= 0 or k < 1: return 0 if k > n: k = n if n < __q_cache.shape[0]: return __q_cache[n][k] return log_q_approx(n, k)
[docs]def get_v(u, epsilon=1e-8): """get_v Parameters ---------- u : int epsilon : float Returns ------- """ v = u delta = 1 while delta > epsilon: n_v = u * np.sqrt(spence(np.exp(-v))) delta = abs(n_v - v) v = n_v return v
[docs]def log_q_approx_small(n, k): """log_q_approx_small Parameters ---------- n : int k : int Returns ------- """ return lbinom(n - 1, k - 1) - loggamma(k + 1)
[docs]def log_q_approx(n, k): """log_q_approx Parameters ---------- n : int k : int Returns ------- """ if k < pow(n, 1 / 4.): return log_q_approx_small(n, k) u = k / np.sqrt(n) v = get_v(u) lf = np.log(v) - np.log1p(- np.exp(-v) * (1 + u * u / 2)) / 2 - np.log(2) * 3 / 2. - np.log(u) - np.log(np.pi) g = 2 * v / u - u * np.log1p(-np.exp(-v)) return lf - np.log(n) + np.sqrt(n) * g
[docs]def init_q_cache(n_max, __q_cache=np.array([], ndmin=2)): """Initiate the look-up table for :math:q(m, n). Parameters ---------- n_max : int __q_cache : :class:numpy.ndarray (required, default: np.array([], ndmin=2)) Returns ------- """ old_n = __q_cache.shape[0] if old_n >= n_max: return __q_cache = np.resize(__q_cache, [n_max + 1, n_max + 1]) __q_cache.fill(-np.inf) return __fill_cache(__q_cache, n_max)
@njit(cache=True) def __fill_cache(cache, n_max): for n in range(1, n_max + 1): cache[n][1] = 0 for k in range(2, n + 1): cache[n][k] = log_sum(cache[n][k], cache[n][k - 1]) if n > k: cache[n][k] = log_sum(cache[n][k], cache[n - k][k]) return cache
[docs]@njit(cache=True) def log_sum(a, b): """log_sum Parameters ---------- a : int b : int Returns ------- """ return np.maximum(a, b) + np.log1p(np.exp(-np.abs(a - b)))
[docs]def lbinom(n, k): """Return log of binom(n, k).""" if type(n) in [float, int, np.int64, np.float64]: n = np.array([n]) k = np.array([k]) return (gammaln(np.array([float(x) for x in n + 1])) - gammaln(np.array([float(x) for x in n - k + 1])) - gammaln(np.array([float(x) for x in k + 1])))