Implement UncertainBaseState.collect_marginal_multigraph() and marginal_multigraph_entropy()

src/graph/inference/uncertain/graph_blockmodel_uncertain.cc

@@ -63,6 +63,12 @@ void collect_xmarginal_dispatch(GraphInterface& gi, GraphInterface& ui,
                                 boost::any aecount, boost::any ax,
                                 boost::any axsum, boost::any ax2sum);
 
+void collect_marginal_count_dispatch(GraphInterface& gi, GraphInterface& ui,
+                                     boost::any aecount, boost::any aexs,
+                                     boost::any aexc);
+
+void marginal_count_entropy(GraphInterface& gi, boost::any aexc, boost::any aeh);
+
 void export_uncertain_state()
 {
     using namespace boost::python;
@@ -115,4 +121,6 @@ void export_uncertain_state()
 
     def("collect_marginal", &collect_marginal_dispatch);
     def("collect_xmarginal", &collect_xmarginal_dispatch);
+    def("collect_marginal_count", &collect_marginal_count_dispatch);
+    def("marginal_count_entropy", &marginal_count_entropy);
 }

src/graph/inference/uncertain/graph_blockmodel_uncertain_marginal.cc

@@ -21,15 +21,6 @@
 using namespace boost;
 using namespace graph_tool;
 
-class dummy_property
-{};
-
-template <class Key>
-constexpr double get(dummy_property&, Key&&) { return 0.; }
-
-template <class Key>
-constexpr void put(dummy_property&, Key&&, double) {}
-
 void collect_marginal_dispatch(GraphInterface& gi, GraphInterface& ui,
                                boost::any aecount)
 {
@@ -38,6 +29,8 @@ void collect_marginal_dispatch(GraphInterface& gi, GraphInterface& ui,
 
     gt_dispatch<>()
         ([&](auto& g, auto& u) { collect_marginal(g, u, ecount,
+                                                  dummy_property(),
+                                                  dummy_property(),
                                                   dummy_property(),
                                                   dummy_property(),
                                                   dummy_property()); },
@@ -59,7 +52,59 @@ void collect_xmarginal_dispatch(GraphInterface& gi, GraphInterface& ui,
 
     gt_dispatch<>()
         ([&](auto& g, auto& u) { collect_marginal(g, u, ecount,
-                                                  x, xsum, x2sum); },
+                                                  x, xsum, x2sum,
+                                                  dummy_property(),
+                                                  dummy_property()); },
          all_graph_views(), all_graph_views())(gi.get_graph_view(),
                                                ui.get_graph_view());
 }
+
+void collect_marginal_count_dispatch(GraphInterface& gi, GraphInterface& ui,
+                                     boost::any aex, boost::any aexs,
+                                     boost::any aexc)
+{
+    typedef eprop_map_t<int32_t>::type ecmap_t;
+    auto ex = any_cast<ecmap_t>(aex);
+
+    typedef eprop_map_t<std::vector<int32_t>>::type emap_t;
+    auto exs = any_cast<emap_t>(aexs);
+    auto exc = any_cast<emap_t>(aexc);
+
+    gt_dispatch<>()
+        ([&](auto& g, auto& u) { collect_marginal(g, u,
+                                                  dummy_property(),
+                                                  ex,
+                                                  dummy_property(),
+                                                  dummy_property(),
+                                                  exs, exc); },
+         all_graph_views(), all_graph_views())(gi.get_graph_view(),
+                                               ui.get_graph_view());
+}
+
+void marginal_count_entropy(GraphInterface& gi, boost::any aexc, boost::any aeh)
+{
+    typedef eprop_map_t<double>::type ehmap_t;
+    auto eh = any_cast<ehmap_t>(aeh);
+
+    gt_dispatch<>()
+        ([&](auto& g, auto exc)
+         {
+             for (auto e : edges_range(g))
+             {
+                 auto& S = eh[e];
+                 S = 0;
+                 size_t N = 0;
+                 for (auto n : exc[e])
+                 {
+                     S -= xlogx_fast(n);
+                     N += n;
+                 }
+                 if (N == 0)
+                     continue;
+                 S /= N;
+                 S += safelog_fast(N);
+             }
+         },
+        all_graph_views(), edge_scalar_vector_properties())
+        (gi.get_graph_view(), aexc);
+}

src/graph/inference/uncertain/graph_blockmodel_uncertain_marginal.hh

@@ -31,9 +31,20 @@ namespace graph_tool
 using namespace boost;
 using namespace std;
 
-template <class Graph, class UGraph, class Eprop, class Xprop>
-void collect_marginal(Graph& g, UGraph& u, Eprop ecount, Xprop x, Xprop xsum,
-                      Xprop x2sum)
+class dummy_property
+{};
+
+template <class Key>
+constexpr double get(dummy_property&, Key&&) { return 0.; }
+
+template <class Key, class Val>
+constexpr void put(dummy_property&, Key&&, Val&&) {}
+
+
+template <class Graph, class UGraph, class Eprop, class Xprop, class XSprop, class Cprop>
+void collect_marginal(Graph& g, UGraph& u, Eprop ecount, Xprop x, XSprop xsum,
+                      XSprop x2sum, [[maybe_unused]] Cprop xs,
+                      [[maybe_unused]] Cprop xcount)
 {
     typedef typename graph_traits<Graph>::edge_descriptor edge_t;
     typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
@@ -57,7 +68,7 @@ void collect_marginal(Graph& g, UGraph& u, Eprop ecount, Xprop x, Xprop xsum,
         {
             ge = add_edge(get<0>(vs), get<1>(vs), g).first;
             emap[vs] = ge;
-            ecount[ge] = 0;
+            put(ecount, ge, 0);
             put(xsum, ge, 0);
             put(x2sum, ge, 0);
         }
@@ -65,9 +76,22 @@ void collect_marginal(Graph& g, UGraph& u, Eprop ecount, Xprop x, Xprop xsum,
         {
             ge = iter->second;
         }
-        ecount[ge]++;
+        put(ecount, ge, get(ecount, ge) + 1);
         put(xsum, ge, get(xsum, ge) + get(x, e));
         put(x2sum, ge, get(x2sum, ge) + get(x, e) * get(x, e));
+        if constexpr (!std::is_same_v<Cprop, dummy_property>)
+        {
+            auto xe = get(x, e);
+            auto& xs_e = xs[ge];
+            auto& xc_e = xcount[ge];
+            auto iter = std::lower_bound(xs_e.begin(), xs_e.end(), xe);
+            if (iter == xs_e.end() || *iter != xe)
+            {
+                iter = xs_e.insert(iter, xe);
+                xc_e.insert(xc_e.begin() + (iter - xs_e.begin()), 0);
+            }
+            xc_e[iter - xs_e.begin()]++;
+        }
     }
 }
 

src/graph_tool/inference/__init__.py

@@ -78,6 +78,7 @@ Auxiliary functions
    ~graph_tool.inference.blockmodel.mf_entropy
    ~graph_tool.inference.blockmodel.bethe_entropy
    ~graph_tool.inference.blockmodel.microstate_entropy
+   ~graph_tool.inference.blockmodel.marginal_multigraph_entropy
    ~graph_tool.inference.overlap_blockmodel.half_edge_graph
    ~graph_tool.inference.overlap_blockmodel.get_block_edge_gradient
 
@@ -180,6 +181,7 @@ __all__ = ["minimize_blockmodel_dl",
            "mf_entropy",
            "bethe_entropy",
            "microstate_entropy",
+           "marginal_multigraph_entropy",
            "PartitionHist",
            "BlockPairHist",
            "half_edge_graph",

src/graph_tool/inference/uncertain_blockmodel.py

@@ -294,6 +294,49 @@ class UncertainBaseState(object):
         g.ep.eprob.fa /= g.gp.count
         return g
 
+    def collect_marginal_multigraph(self, g=None):
+        r"""Collect marginal latent multigraph during MCMC runs.
+
+        Parameters
+        ----------
+        g : :class:`~graph_tool.Graph` (optional, default: ``None``)
+            Previous marginal multigraph.
+
+        Returns
+        -------
+        g : :class:`~graph_tool.Graph`
+            New marginal graph, with internal edge
+            :class:`~graph_tool.EdgePropertyMap` ``"w"`` and ``"wcount"``,
+            containing the edge multiplicities and their respective counts.
+
+        Notes
+        -----
+
+        The mean posterior marginal multiplicity distribution of a multi-edge
+        :math:`(i,j)` is defined as
+
+        .. math::
+
+           \pi_{ij}(w) = \sum_{\boldsymbol G}\delta_{w,G_{ij}}P(\boldsymbol G|\boldsymbol D)
+
+        where :math:`P(\boldsymbol G|\boldsymbol D)` is the posterior
+        probability of a multigraph :math:`\boldsymbol G` given the data.
+
+        """
+
+        if g is None:
+            g = Graph(directed=self.g.is_directed())
+            g.add_vertex(self.g.num_vertices())
+            g.ep.w = g.new_ep("vector<int>")
+            g.ep.wcount = g.new_ep("vector<int>")
+
+        libinference.collect_marginal_count(g._Graph__graph,
+                                            self.u._Graph__graph,
+                                            _prop("e", self.u, self.eweight),
+                                            _prop("e", g, g.ep.w),
+                                            _prop("e", g, g.ep.wcount))
+        return g
+
 class UncertainBlockState(UncertainBaseState):
     r"""Inference state of an uncertain graph, using the stochastic block model as a
     prior.
@@ -470,65 +513,6 @@ class LatentMultigraphBlockState(UncertainBaseState):
                                                  self._state,
                                                  _get_rng())
 
-    def collect_marginal(self, g=None):
-        r"""Collect marginal inferred network during MCMC runs.
-
-        Parameters
-        ----------
-        g : :class:`~graph_tool.Graph` (optional, default: ``None``)
-            Previous marginal graph.
-
-        Returns
-        -------
-        g : :class:`~graph_tool.Graph`
-            New marginal graph, with internal edge
-            :class:`~graph_tool.EdgePropertyMap` ``"x"`` and ``"xdev"``,
-            containing the marginal mean and standard deviation of edge
-            multiplicities, respectively.
-
-        Notes
-        -----
-        The mean posterior marginal multiplicity of an edge :math:`(i,j)` is
-        defined as
-
-        .. math::
-
-           w_{ij} = \sum_{\boldsymbol A}A_{ij}P(\boldsymbol A|\boldsymbol D)
-
-        and likewise the variance is
-
-        .. math::
-
-           \sigma^2_{ij} = \sum_{\boldsymbol A}(A_{ij}-w_{ij})^2P(\boldsymbol A|\boldsymbol D)
-
-        where :math:`P(\boldsymbol A|\boldsymbol D)` is the posterior
-        probability given the data.
-
-        """
-
-        if g is None:
-            g = Graph(directed=self.g.is_directed())
-            g.add_vertex(self.g.num_vertices())
-            g.gp.count = g.new_gp("int", 0)
-            g.ep.count = g.new_ep("int")
-            g.ep.xsum = g.new_ep("double")
-            g.ep.x2sum = g.new_ep("double")
-            g.ep.x = g.new_ep("double")
-            g.ep.xdev = g.new_ep("double")
-
-        u = self.get_graph()
-        x = self.eweight.copy("double")
-        libinference.collect_xmarginal(g._Graph__graph,
-                                       u._Graph__graph,
-                                       _prop("e", u, x),
-                                       _prop("e", g, g.ep.count),
-                                       _prop("e", g, g.ep.xsum),
-                                       _prop("e", g, g.ep.x2sum))
-        g.gp.count += 1
-        g.ep.x.fa = g.ep.xsum.fa / g.gp.count
-        g.ep.xdev.fa = sqrt(g.ep.x2sum.fa / g.gp.count - g.ep.x.fa ** 2)
-        return g
-
 class MeasuredBlockState(UncertainBaseState):
     r"""Inference state of a measured graph, using the stochastic block model as a
     prior.
@@ -1495,3 +1479,47 @@ class PseudoCIsingBlockState(IsingBaseBlockState):
     def _mcmc_sweep_h(self, mcmc_state):
         return libinference.mcmc_pseudo_cising_sweep_h(mcmc_state, self._state,
                                                        _get_rng())
+
+
+def marginal_multigraph_entropy(g, ecount):
+    r"""Compute the entropy of the marginal latent multigraph distribution.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Marginal multigraph.
+    ecount : :class:`~graph_tool.EdgePropertyMap`
+        Vector-valued edge property map containing the counts of edge
+        multiplicities.
+
+    Returns
+    -------
+    eh : :class:`~graph_tool.EdgePropertyMap`
+        Marginal entropy of edge multiplicities.
+
+    Notes
+    -----
+
+    The mean posterior marginal multiplicity distribution of a multi-edge
+    :math:`(i,j)` is defined as
+
+    .. math::
+
+       \pi_{ij}(w) = \sum_{\boldsymbol G}\delta_{w,G_{ij}}P(\boldsymbol G|\boldsymbol D)
+
+    where :math:`P(\boldsymbol G|\boldsymbol D)` is the posterior
+    probability of a multigraph :math:`\boldsymbol G` given the data.
+
+    The corresponding entropy is therefore given (in nats) by
+
+    .. math::
+
+       \mathcal{S}_{ij} = -\sum_w\pi_{ij}(w)\ln \pi_{ij}(w).
+
+    """
+
+    eh = g.new_ep("double")
+    libinference.marginal_count_entropy(g._Graph__graph,
+                                        _prop("e", g, ecount),
+                                        _prop("e", g, eh))
+    return eh