Implement VICenterState

src/graph/inference/Makefile.am

@@ -28,6 +28,9 @@ libgraph_tool_inference_la_SOURCES = \
     blockmodel/graph_blockmodel_multicanonical.cc \
     blockmodel/graph_blockmodel_multicanonical_multiflip.cc \
     blockmodel/graph_blockmodel_multiflip_mcmc.cc \
+    partition_centroid/graph_partition_centroid.cc \
+    partition_centroid/graph_partition_centroid_mcmc.cc \
+    partition_centroid/graph_partition_centroid_multiflip_mcmc.cc \
     overlap/graph_blockmodel_overlap.cc \
     overlap/graph_blockmodel_overlap_exhaustive.cc \
     overlap/graph_blockmodel_overlap_gibbs.cc \
@@ -96,6 +99,7 @@ libgraph_tool_inference_la_include_HEADERS = \
     blockmodel/graph_blockmodel_partition.hh \
     blockmodel/graph_blockmodel_util.hh \
     blockmodel/graph_blockmodel_weights.hh \
+    partition_centroid/graph_partition_centroid.hh \
     overlap/graph_blockmodel_overlap.hh \
     overlap/graph_blockmodel_overlap_mcmc_bundled.hh \
     overlap/graph_blockmodel_overlap_util.hh \

src/graph/inference/graph_inference.cc

@@ -124,6 +124,9 @@ extern void export_pseudo_cising_mcmc_h();
 extern void export_pseudo_ising_state();
 extern void export_pseudo_ising_mcmc();
 extern void export_pseudo_ising_mcmc_h();
+extern void export_vi_center_state();
+extern void export_vi_center_mcmc();
+extern void export_vi_multiflip_mcmc();
 
 BOOST_PYTHON_MODULE(libgraph_tool_inference)
 {
@@ -184,6 +187,9 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
     export_pseudo_ising_state();
     export_pseudo_ising_mcmc();
     export_pseudo_ising_mcmc_h();
+    export_vi_center_state();
+    export_vi_center_mcmc();
+    export_vi_multiflip_mcmc();
 
     def("vector_map", vector_map<int32_t>);
     def("vector_map64", vector_map<int64_t>);

src/graph/inference/partition_centroid/graph_partition_centroid.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_tool.hh"
+#include "random.hh"
+
+#include <boost/python.hpp>
+
+#include "graph_partition_centroid.hh"
+#include "../support/graph_state.hh"
+
+using namespace boost;
+using namespace graph_tool;
+
+GEN_DISPATCH(block_state, VICenterState, BLOCK_STATE_params)
+
+python::object make_vi_center_state(boost::python::object ostate)
+{
+    python::object state;
+    block_state::make_dispatch(ostate,
+                               [&](auto& s){state = python::object(s);});
+    return state;
+}
+
+void export_vi_center_state()
+{
+    using namespace boost::python;
+    def("make_vi_center_state", &make_vi_center_state);
+
+    block_state::dispatch
+        ([&](auto* s)
+         {
+             typedef typename std::remove_reference<decltype(*s)>::type state_t;
+             void (state_t::*move_vertex)(size_t, size_t) =
+                 &state_t::move_vertex;
+             double (state_t::*virtual_move)(size_t, size_t, size_t) =
+                 &state_t::virtual_move;
+
+             class_<state_t>
+                 c(name_demangle(typeid(state_t).name()).c_str(),
+                   no_init);
+             c.def("move_vertex", move_vertex)
+                 .def("virtual_move", virtual_move)
+                 .def("entropy", &state_t::entropy);
+         });
+}

src/graph/inference/partition_centroid/graph_partition_centroid.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_PARTITION_CENTROID_HH
+#define GRAPH_PARTITION_CENTROID_HH
+
+#include "config.h"
+
+#include <vector>
+
+#include "../blockmodel/graph_blockmodel_util.hh"
+#include "../support/graph_state.hh"
+
+#include "openmp_lock.hh"
+
+namespace graph_tool
+{
+using namespace boost;
+using namespace std;
+
+typedef multi_array_ref<int32_t,2> bs_t;
+typedef multi_array_ref<int32_t,1> b_t;
+
+#define BLOCK_STATE_params                                                     \
+    ((g, &, all_graph_views, 1))                                               \
+    ((_abg, &, boost::any&, 0))                                                \
+    ((bs,, bs_t, 0))                                                           \
+    ((b,, b_t, 0))
+
+GEN_STATE_BASE(VICenterStateBase, BLOCK_STATE_params)
+
+template <class... Ts>
+class VICenterState
+    : public VICenterStateBase<Ts...>
+{
+public:
+    GET_PARAMS_USING(VICenterStateBase<Ts...>, BLOCK_STATE_params)
+    GET_PARAMS_TYPEDEF(Ts, BLOCK_STATE_params)
+
+    template <class... ATs,
+              typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
+    VICenterState(ATs&&... args)
+        : VICenterStateBase<Ts...>(std::forward<ATs>(args)...),
+        _bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
+        _mrs(_bs.shape()[0]),
+        _nr(_bs.shape()[0]),
+        _N(_bs.shape()[1]),
+        _wr(_N),
+        _empty_pos(_N),
+        _candidate_pos(_N),
+        _bclabel(_N),
+        _pclabel(_N)
+    {
+        for (size_t r : _b)
+            _wr[r]++;
+
+        for (size_t r = 0; r < _N; ++r)
+        {
+            if (_wr[r] == 0)
+                add_element(_empty_blocks, _empty_pos, r);
+            else
+                add_element(_candidate_blocks, _candidate_pos, r);
+        }
+
+        for (size_t i = 0; i < _mrs.size(); ++i)
+        {
+            for (size_t v = 0; v < _N; ++v)
+            {
+                auto r = _b[v];
+                auto s = _bs[i][v];
+                _mrs[i][{r,s}]++;
+                _nr[i][s]++;
+            }
+        }
+    }
+
+    typedef typename
+        std::conditional<is_directed_::apply<g_t>::type::value,
+                         GraphInterface::multigraph_t,
+                         undirected_adaptor<GraphInterface::multigraph_t>>::type
+        bg_t;
+    bg_t& _bg;
+
+    std::vector<gt_hash_map<std::tuple<size_t, size_t>, size_t>> _mrs;
+    std::vector<gt_hash_map<size_t, size_t>> _nr;
+
+    size_t _N;
+
+    std::vector<size_t> _wr;
+
+    std::vector<size_t> _empty_blocks;
+    std::vector<size_t> _empty_pos;
+    std::vector<size_t> _candidate_blocks;
+    std::vector<size_t> _candidate_pos;
+
+    std::vector<size_t> _bclabel;
+    std::vector<size_t> _pclabel;
+
+    constexpr static BlockStateVirtualBase* _coupled_state = nullptr;
+
+    typedef int m_entries_t;
+
+    bool _egroups_update = true;
+
+    // =========================================================================
+    // State modification
+    // =========================================================================
+
+    void move_vertex(size_t v, size_t nr)
+    {
+        size_t r = _b[v];
+        if (nr == r)
+            return;
+
+        _wr[r]--;
+        _wr[nr]++;
+
+        for (size_t i = 0; i < _mrs.size(); ++i)
+        {
+            auto& mrsi = _mrs[i];
+            size_t s = _bs[i][v];
+            auto iter = mrsi.find({r, s});
+            assert(iter != mrsi.end());
+            iter->second--;
+            if (iter->second == 0)
+                mrsi.erase(iter);
+            mrsi[{nr,s}]++;
+        }
+
+        if (_wr[r] == 0)
+        {
+            add_element(_empty_blocks, _empty_pos, r);
+            remove_element(_candidate_blocks, _candidate_pos, r);
+        }
+
+        if (_wr[nr] == 1)
+        {
+            remove_element(_empty_blocks, _empty_pos, nr);
+            add_element(_candidate_blocks, _candidate_pos, nr);
+        }
+
+        _b[v] = nr;
+    }
+
+    size_t virtual_remove_size(size_t v)
+    {
+        return _wr[_b[v]] - 1;
+    }
+
+    constexpr size_t add_block()
+    {
+        return 0;
+    }
+
+    double virtual_move(size_t v, size_t r, size_t nr)
+    {
+        if (r == nr)
+            return 0;
+
+        double Sb = 0;
+        double Sa = 0;
+
+        for (size_t i = 0; i < _mrs.size(); ++i)
+        {
+            auto& mrsi = _mrs[i];
+            size_t s = _bs[i][v];
+
+            size_t mrs = mrsi[{r,s}];
+            assert(mrs > 0);
+            auto iter = mrsi.find({nr, s});
+            size_t mnrs = (iter != mrsi.end()) ? iter->second : 0;
+
+            Sb += -2 * (xlogx_fast(mrs) + xlogx_fast(mnrs));
+            Sa += -2 * (xlogx_fast(mrs-1) + xlogx_fast(mnrs+1));
+
+            Sb += xlogx_fast(_wr[r]) + xlogx_fast(_wr[nr]);
+            Sa += xlogx_fast(_wr[r]-1) + xlogx_fast(_wr[nr]+1);
+        }
+        return (Sa - Sb);
+    }
+
+    size_t get_empty_block(size_t)
+    {
+        return _empty_blocks.back();
+    }
+
+    size_t sample_block(size_t, double, double d, rng_t& rng)
+    {
+        std::bernoulli_distribution new_r(d);
+        if (d > 0 && !_empty_blocks.empty() && new_r(rng))
+            return uniform_sample(_empty_blocks, rng);
+        return uniform_sample(_candidate_blocks, rng);
+    }
+
+    // Computes the move proposal probability
+    double get_move_prob(size_t, size_t r, size_t s, double, double d, bool reverse)
+    {
+        size_t B = _candidate_blocks.size();
+        if (reverse)
+        {
+            if (_wr[s] == 1)
+                return d;
+            if (_wr[r] == 0)
+                B++;
+        }
+        else
+        {
+            if (_wr[s] == 0)
+                return d;
+        }
+
+        if (B == _N)
+            d = 0;
+
+        return (1. - d) / B;
+    }
+
+    template <class MEntries>
+    double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
+                         bool reverse, MEntries&&)
+    {
+        return get_move_prob(v, r, s, c, d, reverse);
+    }
+
+    template <class EArgs>
+    double virtual_move(size_t v, size_t r, size_t nr, EArgs&&)
+    {
+        return virtual_move(v, r, nr);
+    }
+
+    template <class EArgs, class MEntries>
+    double virtual_move(size_t v, size_t r, size_t nr, EArgs&&, MEntries&&)
+    {
+        return virtual_move(v, r, nr);
+    }
+
+    double entropy()
+    {
+        double S = 0, S_n = 0;
+        for (auto nr : _wr)
+            S_n += xlogx_fast(nr);
+        for (size_t i = 0; i < _mrs.size(); ++i)
+        {
+            for (auto c : _mrs[i])
+                S -= 2 * xlogx_fast(c.second);
+            for (auto rn : _nr[i])
+                S += xlogx_fast(rn.second);
+            S += S_n;
+        }
+        return S;
+    }
+
+    void init_mcmc(double, double)
+    {
+    }
+
+    constexpr size_t node_weight(size_t)
+    {
+        return 1;
+    }
+
+    bool is_last(size_t v)
+    {
+        return _wr[_b[v]] == 1;
+    }
+
+    constexpr bool allow_move(size_t, size_t)
+    {
+        return true;
+    }
+
+};
+
+} // graph_tool namespace
+
+#endif //GRAPH_PARTITION_CENTROID_HH

src/graph/inference/partition_centroid/graph_partition_centroid_mcmc.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_tool.hh"
+#include "random.hh"
+
+#include <boost/python.hpp>
+
+#include "graph_partition_centroid.hh"
+#include "../blockmodel/graph_blockmodel_mcmc.hh"
+#include "../loops/mcmc_loop.hh"
+
+using namespace boost;
+using namespace graph_tool;
+
+GEN_DISPATCH(block_state, VICenterState, BLOCK_STATE_params)
+
+template <class State>
+GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
+             MCMC_BLOCK_STATE_params(State))
+
+python::object vi_mcmc_sweep(python::object omcmc_state,
+                             python::object oblock_state,
+                             rng_t& rng)
+{
+    python::object ret;
+    auto dispatch = [&](auto& block_state)
+    {
+        typedef typename std::remove_reference<decltype(block_state)>::type
+            state_t;
+
+        mcmc_block_state<state_t>::make_dispatch
+           (omcmc_state,
+            [&](auto& s)
+            {
+                auto ret_ = mcmc_sweep(s, rng);
+                ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
+            });
+    };
+    block_state::dispatch(oblock_state, dispatch);
+    return ret;
+}
+
+void export_vi_center_mcmc()
+{
+    using namespace boost::python;
+    def("vi_mcmc_sweep", &vi_mcmc_sweep);
+}

src/graph/inference/partition_centroid/graph_partition_centroid_multiflip_mcmc.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_tool.hh"
+#include "random.hh"
+
+#include <boost/python.hpp>
+
+#include "graph_partition_centroid.hh"
+#include "../blockmodel/graph_blockmodel_multiflip_mcmc.hh"
+#include "../loops/mcmc_loop.hh"
+
+using namespace boost;
+using namespace graph_tool;
+
+GEN_DISPATCH(block_state, VICenterState, BLOCK_STATE_params)
+
+template <class State>
+GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
+             MCMC_BLOCK_STATE_params(State))
+
+python::object vi_multiflip_mcmc_sweep(python::object omcmc_state,
+                                       python::object oblock_state,
+                                       rng_t& rng)
+{
+    python::object ret;
+    auto dispatch = [&](auto& block_state)
+    {
+        typedef typename std::remove_reference<decltype(block_state)>::type
+            state_t;
+
+        mcmc_block_state<state_t>::make_dispatch
+           (omcmc_state,
+            [&](auto& s)
+            {
+                auto ret_ = mcmc_sweep(s, rng);
+                ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
+            });
+    };
+    block_state::dispatch(oblock_state, dispatch);
+    return ret;
+}
+
+void export_vi_multiflip_mcmc()
+{
+    using namespace boost::python;
+    def("vi_multiflip_mcmc_sweep", &vi_multiflip_mcmc_sweep);
+}

src/graph_tool/Makefile.am

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

src/graph_tool/inference/__init__.py

@@ -155,6 +155,7 @@ __all__ = ["minimize_blockmodel_dl",
            "OverlapBlockState",
            "LayeredBlockState",
            "NestedBlockState",
+           "VICenterState",
            "LatentMultigraphBlockState",
            "UncertainBlockState",
            "MeasuredBlockState",
@@ -199,3 +200,4 @@ from . blockmodel_em import *
 from . util import *
 from . modularity import *
 from . latent_multigraph import *
+from . partition_centroid import *

src/graph_tool/inference/partition_centroid.py

+#! /usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# graph_tool -- a general graph manipulation python module
+#
+# 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/>.
+
+from __future__ import division, absolute_import, print_function
+import sys
+if sys.version_info < (3,):
+    range = xrange
+
+from .. import _degree, _prop, Graph, GraphView, _get_rng, Vector_size_t
+from . blockmodel import DictState, get_entropy_args, _bm_test
+
+from .. dl_import import dl_import
+dl_import("from . import libgraph_tool_inference as libinference")
+
+import numpy as np
+import math
+
+class VICenterState(object):
+    r"""Obtain the center of a set of partitions, according to the variation of
+    information metric.
+
+    Parameters
+    ----------
+    bs : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
+        Initial block labels on the vertices. If not supplied, it will be
+        randomly sampled.
+    B : ``int`` (optional, default: ``None``)
+        Number of blocks (or vertex groups). If not supplied it will be obtained
+        from the parameter ``b``.
+    """
+
+    def __init__(self, bs, b=None):
+
+        self.bs = bs = np.asarray(bs, dtype="int32")