Implement generate_maxent_sbm()

src/graph/generation/Makefile.am

@@ -25,6 +25,7 @@ libgraph_tool_generation_la_SOURCES = \
     graph_geometric.cc \
     graph_lattice.cc \
     graph_line_graph.cc \
+    graph_maxent_sbm.cc \
     graph_predecessor.cc \
     graph_price.cc \
     graph_rewiring.cc \
@@ -42,6 +43,7 @@ libgraph_tool_generation_la_include_HEADERS = \
     graph_generation.hh \
     graph_geometric.hh \
     graph_lattice.hh \
+    graph_maxent_sbm.hh \
     graph_predecessor.hh \
     graph_price.hh \
     graph_rewiring.hh \

src/graph/generation/graph_generation.cc

@@ -122,6 +122,8 @@ void community_network_eavg(GraphInterface& gi, GraphInterface& cgi,
                             boost::any eweight, boost::python::list aeprops,
                             bool self_loops, bool parallel_edges);
 
+void export_maxent_sbm();
+
 using namespace boost::python;
 
 BOOST_PYTHON_MODULE(libgraph_tool_generation)
@@ -144,6 +146,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_generation)
     def("community_network", &community_network);
     def("community_network_vavg", &community_network_vavg);
     def("community_network_eavg", &community_network_eavg);
+    export_maxent_sbm();
 
     class_<Sampler<int, boost::mpl::false_>>("Sampler",
                                              init<const vector<int>&, const vector<double>&>())

src/graph/generation/graph_maxent_sbm.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2019 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_maxent_sbm.hh"
+#include "numpy_bind.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void generate_maxent_sbm(GraphInterface& gi, boost::any ab,
+                         boost::python::object ors, boost::python::object oss,
+                         boost::python::object omrs, boost::any ain_theta,
+                         boost::any aout_theta, bool multi, bool self_loops,
+                         rng_t& rng)
+{
+    auto rs = get_array<int64_t, 1>(ors);
+    auto ss = get_array<int64_t, 1>(oss);
+    auto mrs = get_array<double, 1>(omrs);
+
+    typedef vprop_map_t<int32_t>::type bmap_t;
+    auto b = any_cast<bmap_t>(ab).get_unchecked();
+
+    typedef vprop_map_t<double>::type dmap_t;
+    auto in_theta = any_cast<dmap_t>(ain_theta).get_unchecked();
+    auto out_theta = any_cast<dmap_t>(aout_theta).get_unchecked();
+
+    if (multi)
+    {
+        run_action<>()
+            (gi, [&](auto& g) { gen_maxent_sbm<true>(g, b, rs, ss, mrs, in_theta,
+                                                     out_theta, self_loops, rng); })();
+    }
+    else
+    {
+        run_action<>()
+            (gi, [&](auto& g) { gen_maxent_sbm<false>(g, b, rs, ss, mrs, in_theta,
+                                                      out_theta, self_loops, rng); })();
+    }
+}
+
+SBMFugacities get_sbm_fugacities(boost::python::object ors,
+                                 boost::python::object oss,
+                                 boost::python::object oers,
+                                 boost::python::object odegs_in,
+                                 boost::python::object odegs_out,
+                                 boost::python::object ob, bool directed,
+                                 bool multigraph, bool self_loops)
+{
+    auto rs = get_array<int64_t,1>(ors);
+    auto ss = get_array<int64_t,1>(oss);
+    auto ers = get_array<double,1>(oers);
+    auto degs_in = get_array<double,1>(odegs_in);
+    auto degs_out = get_array<double,1>(odegs_out);
+    auto b = get_array<int32_t,1>(ob);
+    return SBMFugacities(rs, ss, ers, degs_in, degs_out, b, directed,
+                         multigraph, self_loops);
+}
+
+using namespace boost::python;
+
+void export_maxent_sbm()
+{
+    def("get_sbm_fugacities", &get_sbm_fugacities);
+    def("gen_maxent_sbm", &generate_maxent_sbm);
+
+    class_<SBMFugacities>("SBMFugacities", no_init)
+        .def("pack", &SBMFugacities::pack)
+        .def("unpack", &SBMFugacities::unpack)
+        .def("get_f", &SBMFugacities::get_f)
+        .def("get_diff", &SBMFugacities::get_diff)
+        .def("new_x", &SBMFugacities::new_x)
+        .def("norm", &SBMFugacities::norm)
+        .def("export_args",
+             +[](SBMFugacities& state, python::object ors, python::object oss,
+                 python::object omrs, python::object odegs_in,
+                 python::object odegs_out, python::object otheta_in,
+                 python::object otheta_out, python::object ob)
+              {
+                  auto rs = get_array<int64_t,1>(ors);
+                  auto ss = get_array<int64_t,1>(oss);
+                  auto mrs = get_array<double,1>(omrs);
+                  auto degs_in = get_array<double,1>(odegs_in);
+                  auto degs_out = get_array<double,1>(odegs_out);
+                  auto theta_in = get_array<double,1>(otheta_in);
+                  auto theta_out = get_array<double,1>(otheta_out);
+                  auto b = get_array<int32_t,1>(ob);
+                  state.export_args(rs, ss, mrs, degs_in, degs_out, theta_in,
+                                    theta_out, b);
+              });
+}

src/graph/generation/graph_maxent_sbm.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2019 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_MAXENT_SBM_HH
+#define GRAPH_MAXENT_SBM_HH
+
+#include <tuple>
+#include <iostream>
+#include <boost/functional/hash.hpp>
+#include <boost/multi_array.hpp>
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+#include "sampler.hh"
+#include "urn_sampler.hh"
+
+#include "random.hh"
+
+#include "hash_map_wrap.hh"
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+
+class SBMFugacities
+{
+public:
+
+    template <class IVec, class FVec, class Vec, class Bvec>
+    SBMFugacities(IVec&& rs, IVec&& ss, FVec&& ers, Vec&& degs_in, Vec&& degs_out,
+                  Bvec&& b, bool directed, bool multigraph, bool self_loops)
+        : _directed(directed), _multigraph(multigraph), _self_loops(self_loops)
+    {
+
+        size_t N = degs_in.size();
+        double E2 = 0;
+        _B = *std::max_element(b.begin(), b.end()) + 1;
+        for (size_t i = 0; i < rs.size(); ++i)
+        {
+            size_t r = rs[i];
+            size_t s = ss[i];
+            if (r >= _mrs.size())
+            {
+                _mrs.resize(r + 1);
+                _ers.resize(r + 1);
+            }
+            if (s >= _msr.size())
+                _msr.resize(s + 1);
+            E2 += _mrs[r][s] = _msr[s][r] = _ers[r][s] = ers[i];
+            _B = std::max(_B, std::max(r, s) + 1);
+        }
+
+        _mrs.resize(_B);
+        _msr.resize(_B);
+        _ers.resize(_B);
+
+        if (multigraph)
+        {
+            for (auto& mr : _mrs)
+                for (auto& m : mr)
+                    m.second /= E2;
+            for (auto& mr : _msr)
+                for (auto& m : mr)
+                    m.second /= E2;
+        }
+
+        std::vector<gt_hash_map<double, size_t>> vertices_in(_B),
+            vertices_out(_B);
+        std::vector<size_t> nr(_B);
+
+        for (size_t i = 0; i < N; ++i)
+        {
+            vertices_in[b[i]][degs_in[i]]++;
+            vertices_out[b[i]][degs_out[i]]++;
+            nr[b[i]]++;
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            _rtheta_in.emplace_back();
+            _rdegs_in.emplace_back();
+            _in_pos.emplace_back();
+
+            double S = 0;
+            for (auto& rc : vertices_in[r])
+            {
+                double t = rc.first;
+                _in_pos[r][t] = _rtheta_in[r].size();
+                if (multigraph && t > sqrt(nr[r]-1))
+                    t = sqrt(nr[r]-1);
+                _rtheta_in[r].emplace_back(t, rc.second);
+                _rdegs_in[r].push_back(rc.first);
+                S += t * rc.second;
+            }
+            if (multigraph)
+            {
+                for (auto& rc : _rtheta_in[r])
+                    rc.first /= sqrt(S);
+            }
+            else
+            {
+                for (auto& rc : _rtheta_in[r])
+                    rc.first /= S;
+            }
+
+            _rtheta_out.emplace_back();
+            _rdegs_out.emplace_back();
+            _out_pos.emplace_back();
+            S = 0;
+            for (auto rc : vertices_out[r])
+            {
+                double t = rc.first;
+                _out_pos[r][t] = _rtheta_out[r].size();
+                if (multigraph && t > sqrt(nr[r]-1))
+                    t = sqrt(nr[r]-1);
+                _rtheta_out[r].emplace_back(t, rc.second);
+                _rdegs_out[r].push_back(rc.first);
+                S += t * rc.second;
+            }
+            if (multigraph)
+            {
+                for (auto& rc : _rtheta_out[r])
+                    rc.first /= sqrt(S);
+            }
+            else
+            {
+                for (auto& rc : _rtheta_out[r])
+                    rc.first /= sqrt(S);
+            }
+        }
+    }
+
+    void norm()
+    {
+        std::vector<double> t_in(_B), t_out(_B);
+        for (size_t r = 0; r < _B; ++r)
+        {
+            t_in[r] = 0;
+            for (auto& rc : _rtheta_in[r])
+                t_in[r] += rc.first * rc.second;
+            for (auto& rc : _rtheta_in[r])
+                rc.first /= t_in[r];
+            t_out[r] = 0;
+            for (auto& rc : _rtheta_out[r])
+                t_out[r] += rc.first * rc.second;
+            for (auto& rc : _rtheta_out[r])
+                rc.first /= t_out[r];
+        }
+    }
+
+    void pack(std::vector<double>& x)
+    {
+        x.clear();
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& rc : _rtheta_out[r])
+                x.push_back(rc.first);
+            if (_directed)
+            {
+                for (auto& rc : _rtheta_in[r])
+                    x.push_back(rc.first);
+            }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& m : _mrs[r])
+            {
+                auto s = m.first;
+                if (!_directed && r > s)
+                    continue;
+                x.push_back(m.second);
+            }
+        }
+    }
+
+    void unpack(std::vector<double>& x)
+    {
+        size_t pos = 0;
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& rc : _rtheta_out[r])
+                rc.first = x[pos++];
+            if (_directed)
+            {
+                for (auto& rc : _rtheta_in[r])
+                    rc.first = x[pos++];
+            }
+            else
+            {
+                _rtheta_in[r] = _rtheta_out[r];
+            }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+            for (auto& m : _mrs[r])
+            {
+                auto s = m.first;
+                if (!_directed && r > s)
+                    continue;
+                m.second = x[pos++];
+            }
+
+        if (!_directed)
+        {
+            for (size_t r = 0; r < _B; ++r)
+                for (auto& m : _mrs[r])
+                {
+                    auto s = m.first;
+                    if (r > s)
+                        m.second = _mrs[s][r];
+                }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+            for (auto& m : _mrs[r])
+                _msr[m.first][r] = m.second;
+    }
+
+    double get_f()
+    {
+        double L = 0;
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            if (_directed)
+            {
+                for (size_t i = 0; i < _rtheta_in[r].size(); ++i)
+                    L += _rdegs_in[r][i] * log(_rtheta_in[r][i].first) * double(_rtheta_in[r][i].second);
+            }
+            for (size_t i = 0; i < _rtheta_out[r].size(); ++i)
+                L += _rdegs_out[r][i] * log(_rtheta_out[r][i].first) * double(_rtheta_out[r][i].second);
+
+            for (auto& m : _mrs[r])
+            {
+                auto s = m.first;
+                L += _ers[r][s] * log(m.second) / (_directed ? 1. : 2.);
+
+                for (size_t i = 0; i < _rtheta_out[r].size(); ++i)
+                {
+                    double tout = _rtheta_out[r][i].first;
+                    double nout = _rtheta_out[r][i].second;
+
+                    for (auto& tc : _rtheta_in[s])
+                    {
+                        double n = nout * ((s != r || tout != tc.first) ? tc.second : tc.second - 1);
+                        if (!_directed)
+                            n /= 2;
+                        double p = tout * tc.first * _mrs[r][s];
+                        if (!_multigraph)
+                            L -= log1p(p) * n;
+                        else
+                            L += log1p(-p) * n;
+                    }
+                }
+            }
+        }
+
+        return L;
+    }
+
+    void get_diff(std::vector<double>& x)
+    {
+        std::vector<std::vector<std::pair<double, size_t>>> diff_rtheta_in(_rtheta_in),
+            diff_rtheta_out(_rtheta_out);
+
+        vector<gt_hash_map<size_t, double>> diff_mrs(_B);
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (size_t i = 0; i < _rtheta_in[r].size(); ++i)
+            {
+                double tin = _rtheta_in[r][i].first;
+                double S = 0;
+                for (auto& m : _msr[r])
+                {
+                    auto s = m.first;
+                    for (auto& tc : _rtheta_out[s])
+                    {
+                        double n = (s == r && tin == tc.first && !_self_loops)
+                            ? tc.second - 1 : tc.second;
+                        double p;
+                        if (!_multigraph)
+                            p = (tc.first * m.second * n) / (1. + tin * tc.first * m.second);
+                        else
+                            p = (tc.first * m.second * n) / (1. - tin * tc.first * m.second);
+                        S += p;
+                    }
+                }
+                diff_rtheta_in[r][i].first =  _rdegs_in[r][i] / tin - S;
+            }
+
+            for (size_t i = 0; i < _rtheta_out[r].size(); ++i)
+            {
+                double tout = _rtheta_out[r][i].first;
+                double S = 0;
+                for (auto& m : _mrs[r])
+                {
+                    auto s = m.first;
+                    for (auto& tc : _rtheta_in[s])
+                    {
+                        double n = (s == r && tout == tc.first && !_self_loops)
+                            ? tc.second - 1 : tc.second;
+                        double p;
+                        if (!_multigraph)
+                            p = (tc.first * m.second * n) / (1. + tout * tc.first * m.second);
+                        else
+                            p = (tc.first * m.second * n) / (1. - tout * tc.first * m.second);
+                        S += p;
+                    }
+                }
+                diff_rtheta_out[r][i].first = _rdegs_out[r][i] / tout - S;
+            }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& m : _mrs[r])
+            {
+                auto s = m.first;
+                double S = 0;
+                for (auto& tc_out : _rtheta_out[r])
+                    for (auto& tc_in : _rtheta_in[s])
+                    {
+                        double n = (r == s && tc_out.first == tc_in.first
+                                    && !_self_loops) ?
+                            tc_out.second * (tc_in.second - 1) :
+                            tc_out.second * tc_in.second;
+                        double p = tc_out.first * tc_in.first;
+                        if (!_multigraph)
+                            p = p * n / (1. + p * m.second);
+                        else
+                            p = p * n / (1. - p * m.second);
+                        S += p;
+                    }
+                diff_mrs[r][s] = _ers[r][s] / m.second - S;
+                if (r == s)
+                    diff_mrs[r][s] /= 2;
+            }
+        }
+
+        x.clear();
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& rc : diff_rtheta_out[r])
+                x.push_back(rc.first);
+            if (_directed)
+            {
+                for (auto& rc : diff_rtheta_in[r])
+                    x.push_back(rc.first);
+            }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& m : diff_mrs[r])
+            {
+                auto s = m.first;
+                if (!_directed && r > s)
+                    continue;
+                x.push_back(m.second);
+            }
+        }
+    }
+
+    void new_x(std::vector<double>& x)
+    {
+        std::vector<std::vector<std::pair<double, size_t>>> new_rtheta_in(_rtheta_in),
+            new_rtheta_out(_rtheta_out);
+
+        std::vector<gt_hash_map<size_t, double>> new_mrs(_B);
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (size_t i = 0; i < _rtheta_in[r].size(); ++i)
+            {
+                double tin = _rtheta_in[r][i].first;
+                double S = 0;
+                for (auto& m : _msr[r])
+                {
+                    auto s = m.first;
+                    for (auto& tc : _rtheta_out[s])
+                    {
+                        double n = (s == r && tin == tc.first && !_self_loops)
+                            ? tc.second - 1 : tc.second;
+                        if (!_multigraph)
+                        {
+                            S += (tc.first * m.second * n) / (1. + tin * tc.first * m.second);
+                        }
+                        else
+                        {
+                            S += (tc.first * m.second * n) / (1. - tin * tc.first * m.second);
+                        }
+                    }
+                }
+                new_rtheta_in[r][i].first = _rdegs_in[r][i] / S;
+            }
+
+            for (size_t i = 0; i < _rtheta_out[r].size(); ++i)
+            {
+                double tout = _rtheta_out[r][i].first;
+                double S = 0;
+                for (auto& m : _mrs[r])
+                {
+                    auto s = m.first;
+                    for (auto& tc : _rtheta_in[s])
+                    {
+                        double n = (s == r && tout == tc.first && !_self_loops)
+                            ? tc.second - 1 : tc.second;
+                        if (!_multigraph)
+                            S += (tc.first * m.second * n) / (1. + tout * tc.first * m.second);
+                        else
+                            S += (tc.first * m.second * n) / (1. - tout * tc.first * m.second);
+                    }
+                }
+                new_rtheta_out[r][i].first = _rdegs_out[r][i] / S;
+            }
+        }
+
+        for (size_t r = 0; r < _B; ++r)
+        {
+            for (auto& m : _mrs[r])
+            {
+                auto s = m.first;
+                double S = 0;
+                for (auto& tc_out : _rtheta_out[r])
+                    for (auto& tc_in : _rtheta_in[s])
+                    {
+                        double n = (r == s && tc_out.first == tc_in.first
+                                    && !_self_loops) ?
+                            tc_out.second * (tc_in.second - 1) :