inference.blockmodel.BlockState: add 'allow_empty' option

src/graph/inference/graph_blockmodel.hh

@@ -59,7 +59,8 @@ typedef mpl::vector2<simple_degs_t, degs_map_t> degs_tr;
     ((ecov_type,, int, 0))                                                     \
     ((ecov,, eprop_map_t<double>::type, 0))                                    \
     ((bcovsum,, vprop_map_t<double>::type, 0))                                 \
-    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))
+    ((ignore_degrees,, typename vprop_map_t<uint8_t>::type, 0))                \
+    ((allow_empty,, bool, 0))
 
 GEN_STATE_BASE(BlockStateBase, BLOCK_STATE_params)
 
@@ -1214,7 +1215,8 @@ public:
             for (size_t c = 0; c < C; ++c)
                 _partition_stats.emplace_back(_g, _b, vcs[c], E, B,
                                               _vweight, _eweight, _degs,
-                                              _ignore_degrees, _bmap);
+                                              _ignore_degrees, _bmap,
+                                              _allow_empty);
 
             for (auto r : vertices_range(_bg))
                 _partition_stats[rc[r]].get_r(r);

src/graph/inference/graph_blockmodel_util.hh

@@ -149,12 +149,14 @@ public:
               class Mprop, class Vlist>
     partition_stats_t(Graph& g, Vprop& b, Vlist& vlist, size_t E, size_t B,
                       VWprop& vweight, Eprop& eweight, Degs& degs,
-                      const Mprop& ignore_degree, std::vector<size_t>& bmap)
+                      const Mprop& ignore_degree, std::vector<size_t>& bmap,
+                      bool allow_empty)
         : _bmap(bmap)
     {
         _N = 0;
         _E = E;
         _total_B = B;
+        _allow_empty = allow_empty;
 
         for (auto v : vlist)
         {
@@ -212,7 +214,10 @@ public:
     double get_partition_dl()
     {
         double S = 0;
-        S += lbinom(_actual_B + _N - 1, _N);
+        if (_allow_empty)
+            S += lbinom(_actual_B + _N - 1, _N);
+        else
+            S += lbinom(_N - 1, _actual_B - 1);
         S += lgamma(_N + 1);
         for (auto nr : _total)
             S -= lgamma(nr + 1);
@@ -297,8 +302,16 @@ public:
 
         if (dB != 0)
         {
-            S_b += lbinom(_actual_B + _N - 1, _N);
-            S_a += lbinom(_actual_B + dB + _N - 1, _N);
+            if (_allow_empty)
+            {
+                S_b += lbinom(_actual_B + _N - 1, _N);
+                S_a += lbinom(_actual_B + dB + _N - 1, _N);
+            }
+            else
+            {
+                S_b += lbinom(_N - 1, _actual_B - 1);
+                S_a += lbinom(_N - 1, _actual_B + dB - 1);
+            }
         }
 
         return S_a - S_b;
@@ -491,6 +504,7 @@ private:
     size_t _E;
     size_t _actual_B;
     size_t _total_B;
+    bool _allow_empty;
     vector<map_t> _hist;
     vector<int> _total;
     vector<int> _ep;

src/graph_tool/inference/__init__.py

@@ -68,6 +68,7 @@ Auxiliary functions
 .. autosummary::
    :nosignatures:
 
+   model_entropy
    mf_entropy
    bethe_entropy
    half_edge_graph
@@ -119,6 +120,7 @@ __all__ = ["minimize_blockmodel_dl",
            "hierarchy_minimize",
            "EMBlockState",
            "em_infer",
+           "model_entropy",
            "mf_entropy",
            "bethe_entropy",
            "half_edge_graph",

src/graph_tool/inference/blockmodel.py

@@ -100,6 +100,9 @@ class BlockState(object):
     deg_corr : ``bool`` (optional, default: ``True``)
         If ``True``, the degree-corrected version of the blockmodel ensemble will
         be assumed, otherwise the traditional variant will be used.
+    allow_empty : ``bool`` (optional, default: ``True``)
+        If ``True``, partition description length computed will allow for empty
+        groups.
     max_BE : ``int`` (optional, default: ``1000``)
         If the number of blocks exceeds this value, a sparse matrix is used for
         the block graph. Otherwise a dense matrix will be used.
@@ -107,7 +110,7 @@ class BlockState(object):
 
     def __init__(self, g, eweight=None, vweight=None, b=None, B=None, ecov=None,
                  ecov_type=None, clabel=None, pclabel=None, deg_corr=True,
-                 max_BE=1000, **kwargs):
+                 allow_empty=True, max_BE=1000, **kwargs):
         kwargs = kwargs.copy()
 
         # initialize weights to unity, if necessary
@@ -244,6 +247,7 @@ class BlockState(object):
             self.ecov = self.g.new_ep("double")
             self.bcovsum = self.bg.new_vp("double")
 
+        self.allow_empty = allow_empty
         self._abg = self.bg._get_any()
         self._state = libinference.make_block_state(self, _get_rng())
 
@@ -341,8 +345,11 @@ class BlockState(object):
                            vweight=self.wr if vweight else self.bg.new_vp("int", 1),
                            b=self.bg.vertex_index.copy("int") if b is None else b,
                            deg_corr=deg_corr,
+                           allow_empty=kwargs.get("allow_empty",
+                                                  self.allow_empty),
                            degs=degs,
-                           max_BE=self.max_BE, **kwargs)
+                           max_BE=self.max_BE,
+                           **dmask(kwargs, ["allow_empty"]))
         return state
 
     def get_bclabel(self):
@@ -1478,7 +1485,7 @@ class BlockState(object):
                                            "vertex_color",
                                            "edge_gradient"]))
 
-def model_entropy(B, N, E, directed=False, nr=None):
+def model_entropy(B, N, E, directed=False, nr=None, allow_empty=True):
     r"""Computes the amount of information necessary for the parameters of the
     traditional blockmodel ensemble, for ``B`` blocks, ``N`` vertices, ``E``
     edges, and either a directed or undirected graph.
@@ -1505,18 +1512,27 @@ def model_entropy(B, N, E, directed=False, nr=None):
     number of :math:`k` combinations with repetitions from a set of size
     :math:`n`.
 
-    The total information necessary to describe the model is then,
+    For the node partition we assume a two-level Bayesian hierarchy, where first
+    the group size histogram is generated, and conditioned on it the partition,
+    which leads to a description length:
 
     .. math::
 
-       \mathcal{L}_t = \ln\Omega_m + \ln\left(\!\!{B \choose N}\!\!\right) +
-       \ln N! - \sum_r \ln n_r!,
+       \mathcal{L}_p = \ln\left(\!\!{B \choose N}\!\!\right) + \ln N! - \sum_r \ln n_r!,
+
+    where :math:`n_r` is the number of nodes in block :math:`r`. If we forbid
+    empty groups, with ``allow_empty == False``, this changes slightly to
+
+       \mathcal{L}_p = \ln {N - 1 \choose B - 1} + \ln N! - \sum_r \ln n_r!.
 
+    The total information necessary to describe the model is then,
+
+    .. math::
 
-    where the remaining term is the information necessary to describe the block
-    partition, where :math:`n_r` is the number of nodes in block :math:`r`.
+       \mathcal{L}_t = \ln\Omega_m + \mathcal{L}_p.
 
-    If ``nr`` is ``None``, it is assumed :math:`n_r=N/B`.
+    If ``nr`` is ``None``, it is assumed :math:`n_r=N/B`. If ``nr`` is
+    ``False``, the partition term :math:`\mathcal{L}_p` is omitted entirely.
 
     References
     ----------
@@ -1528,6 +1544,10 @@ def model_entropy(B, N, E, directed=False, nr=None):
        structures and high-resolution model selection in large networks ",
        Phys. Rev. X 4, 011047 (2014), :doi:`10.1103/PhysRevX.4.011047`,
        :arxiv:`1310.4377`.
+    .. [peixoto-model-2016] Tiago P. Peixoto, "Model selection and hypothesis
+       testing for large-scale network models with overlapping groups",
+       Phys. Rev. X 5, 011033 (2016), :doi:`10.1103/PhysRevX.5.011033`,
+       :arxiv:`1409.3059`.
 
     """