Implement latent multigraph inference via EM (WIP)

src/graph/inference/Makefile.am

@@ -77,6 +77,7 @@ libgraph_tool_inference_la_SOURCES = \
     support/int_part.cc \
     support/spence.cc \
     graph_inference.cc \
+    graph_latent_multigraph.cc \
     graph_modularity.cc
 
 libgraph_tool_inference_la_include_HEADERS = \
@@ -122,4 +123,5 @@ libgraph_tool_inference_la_include_HEADERS = \
     support/graph_state.hh \
     support/int_part.hh \
     support/util.hh \
-    graph_modularity.hh
+    graph_modularity.hh \
+    graph_latent_multigraph.hh

src/graph/inference/graph_inference.cc

@@ -117,6 +117,7 @@ extern void export_ising_glauber_state();
 extern void export_ising_glauber_mcmc();
 extern void export_marginals();
 extern void export_modularity();
+extern void export_latent_multigraph();
 extern void export_pseudo_cising_state();
 extern void export_pseudo_cising_mcmc();
 extern void export_pseudo_cising_mcmc_h();
@@ -176,6 +177,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
     export_ising_glauber_mcmc();
     export_marginals();
     export_modularity();
+    export_latent_multigraph();
     export_pseudo_cising_state();
     export_pseudo_cising_mcmc();
     export_pseudo_cising_mcmc_h();

src/graph/inference/graph_latent_multigraph.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph_latent_multigraph.hh"
+
+#include <boost/mpl/push_back.hpp>
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void latent_multigraph(GraphInterface& gi, boost::any aw, boost::any atheta_out,
+                       boost::any atheta_in, double epsilon, size_t max_niter)
+{
+    typedef eprop_map_t<double>::type emap_t;
+    typedef vprop_map_t<double>::type vmap_t;
+    auto w = any_cast<emap_t>(aw).get_unchecked();
+    auto theta_out = any_cast<vmap_t>(atheta_out).get_unchecked();
+    auto theta_in = any_cast<vmap_t>(atheta_in).get_unchecked();
+
+    run_action<>()
+        (gi, [&](auto& g){ get_latent_multigraph(g, w, theta_out, theta_in,
+                                                 epsilon, max_niter); })();
+}
+
+using namespace boost::python;
+
+void export_latent_multigraph()
+{
+    def("latent_multigraph", &latent_multigraph);
+}

src/graph/inference/graph_latent_multigraph.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef GRAPH_LATENT_MULTIGRAPH_HH
+#define GRAPH_LATENT_MULTIGRAPH_HH
+
+#include <tuple>
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+
+#include "graph_tool.hh"
+#include "hash_map_wrap.hh"
+
+namespace graph_tool
+{
+
+using namespace std;
+using namespace boost;
+
+template <class Graph, class WMap, class TMap>
+void  get_latent_multigraph(Graph& g, WMap w, TMap theta_out, TMap theta_in,
+                            double epsilon, size_t max_niter)
+{
+    auto wc = w.get_checked();
+    for (auto v : vertices_range(g))
+    {
+        auto tout = theta_out[v];
+        auto tin = theta_in[v];
+        auto e = add_edge(v, v, g).first;
+        if (graph_tool::is_directed(g))
+            wc[e] = tout * tin;
+        else
+            wc[e] = (tout * tin) / 2;
+    }
+
+    double delta = 1 + epsilon;
+    size_t niter = 0;
+    while (delta > epsilon && (niter < max_niter || max_niter == 0))
+    {
+        double M = 0;
+        delta = 0;
+        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH)   \
+            reduction(+:M) reduction(max:delta)
+        parallel_edge_loop_no_spawn
+            (g,
+             [&](auto e)
+             {
+                 auto u = source(e, g);
+                 auto v = target(e, g);
+                 auto l = theta_out[u] * theta_in[v];
+                 auto nw = ((u != v) ? l / (1-exp(-l)) :
+                            (graph_tool::is_directed(g) ? l : l / 2));
+                 auto& ew = w[e];
+                 delta = std::max(delta, abs(nw - ew));
+                 ew = nw;
+                 M += ew;
+             });
+
+        #pragma omp parallel if (num_vertices(g) > OPENMP_MIN_THRESH)   \
+            reduction(max:delta)
+        parallel_vertex_loop_no_spawn
+            (g,
+             [&](auto v)
+             {
+                 if (graph_tool::is_directed(g))
+                 {
+                     auto& tout = theta_out[v];
+                     auto d = out_degreeS()(v, g, w);
+                     auto nt = d / sqrt(M);
+                     delta = std::max(delta, abs(tout - nt));
+                     tout = nt;
+
+                     auto& tin = theta_in[v];
+                     d = in_degreeS()(v, g, w);
+                     nt = d / sqrt(M);
+                     delta = std::max(delta, abs(tin - nt));
+                     tin = nt;
+                 }
+                 else
+                 {
+                     auto& t = theta_out[v];
+                     auto d = out_degreeS()(v, g, w);
+                     auto nt = d / sqrt(2 * M);
+                     delta = std::max(delta, abs(t - nt));
+                     t = nt;
+                 }
+             });
+        niter++;
+    }
+};
+
+} // graph_tool namespace
+
+#endif //GRAPH_LATENT_MULTIGRAPH_HH

src/graph_tool/Makefile.am

@@ -75,6 +75,7 @@ graph_tool_inference_PYTHON = \
     inference/mcmc.py \
     inference/minimize.py \
     inference/modularity.py \
+    inference/latent_multigraph.py \
     inference/util.py
 graph_tool_inferencedir = $(MOD_DIR)/inference
 

src/graph_tool/inference/__init__.py

@@ -175,7 +175,8 @@ __all__ = ["minimize_blockmodel_dl",
            "half_edge_graph",
            "get_block_edge_gradient",
            "get_hierarchy_tree",
-           "modularity"]
+           "modularity",
+           "latent_multigraph"]
 
 from . blockmodel import *
 from . overlap_blockmodel import *
@@ -188,3 +189,4 @@ from . minimize import *
 from . blockmodel_em import *
 from . util import *
 from . modularity import *
+from . latent_multigraph import *

src/graph_tool/inference/latent_multigraph.py

+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# graph_tool -- a general graph manipulation python module
+#
+# Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from .. import _prop
+
+from .. dl_import import dl_import
+dl_import("from . import libgraph_tool_inference as libinference")
+
+from numpy import sqrt
+
+def latent_multigraph(g, epsilon=1e-8, max_niter=10000):
+    r"""
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+
+    Returns
+    -------
+
+    Notes
+    -----
+
+    Examples
+    --------
+    >>> g = gt.collection.data["football"]
+    >>> gt.modularity(g, g.vp.value_tsevans)
+    0.5744393497...
+
+    References
+    ----------
+    """
+
+    g = g.copy()
+    theta_out = g.degree_property_map("out").copy("double")
+    theta_out.fa /= sqrt(theta_out.fa.sum())
+    if g.is_directed():
+        theta_in = g.degree_property_map("in").copy("double")
+        theta_in.fa /= sqrt(theta_in.fa.sum())
+    else:
+        theta_in = theta_out
+
+    w = g.new_ep("double", 1)
+
+    libinference.latent_multigraph(g._Graph__graph,
+                                   _prop("e", g, w),
+                                   _prop("v", g, theta_out),
+                                   _prop("v", g, theta_in),
+                                   epsilon,
+                                   max_niter)
+    return g, w