inference: some improvements to the multicanonical code

src/graph/inference/graph_blockmodel.cc

@@ -139,14 +139,12 @@ boost::python::tuple bethe_entropy(GraphInterface& gi, size_t B, boost::any op,
                                    boost::any opv)
 {
     typedef vprop_map_t<vector<double>>::type vmap_t;
-    typedef eprop_map_t<vector<int32_t>>::type emap_t;
-    emap_t p = any_cast<emap_t>(op);
     vmap_t pv = any_cast<vmap_t>(opv);
 
     double H=0, sH=0, Hmf=0, sHmf=0;
-    run_action<graph_tool::all_graph_views, boost::mpl::true_>()
+    run_action<>()
         (gi,
-         [&](auto& g)
+         [&](auto& g, auto p)
          {
              for (auto v : vertices_range(g))
              {
@@ -204,7 +202,8 @@ boost::python::tuple bethe_entropy(GraphInterface& gi, size_t B, boost::any op,
                      sHmf += pow((log(pi) + 1) * sqrt(pi / sum), 2);
                  }
              }
-         })();
+         },
+         edge_scalar_vector_properties())(op);
 
     return boost::python::make_tuple(H, sH, Hmf, sHmf);
 }

src/graph/inference/graph_blockmodel_layers_multicanonical.cc

@@ -63,7 +63,8 @@ python::object multicanonical_layered_sweep(python::object omulticanonical_state
                       [&](auto& s)
                       {
                           auto ret_ = multicanonical_sweep(s, rng);
-                          ret = python::make_tuple(ret_.first, ret_.second);
+                          ret = python::make_tuple(ret_.first, ret_.second,
+                                                   s._f, s._time, s._refine);
                       });
              },
              false);

src/graph/inference/graph_blockmodel_layers_overlap_multicanonical.cc

@@ -64,7 +64,8 @@ multicanonical_layered_overlap_sweep(python::object omulticanonical_state,
                       [&](auto& s)
                       {
                           auto ret_ = multicanonical_sweep(s, rng);
-                          ret = python::make_tuple(ret_.first, ret_.second);
+                          ret = python::make_tuple(ret_.first, ret_.second,
+                                                   s._f, s._time, s._refine);
                       });
              },
              false);

src/graph/inference/graph_blockmodel_multicanonical.cc

@@ -36,8 +36,8 @@ GEN_DISPATCH(multicanonical_block_state,
              MULTICANONICAL_BLOCK_STATE_params(State))
 
 python::object do_multicanonical_sweep(python::object omulticanonical_state,
-                               python::object oblock_state,
-                               rng_t& rng)
+                                       python::object oblock_state,
+                                       rng_t& rng)
 {
     python::object ret;
     auto dispatch = [&](auto& block_state)
@@ -50,7 +50,8 @@ python::object do_multicanonical_sweep(python::object omulticanonical_state,
             [&](auto& s)
             {
                 auto ret_ = multicanonical_sweep(s, rng);
-                ret = python::make_tuple(ret_.first, ret_.second);
+                ret = python::make_tuple(ret_.first, ret_.second, s._f,
+                                         s._time);
             });
     };
     block_state::dispatch(oblock_state, dispatch);

src/graph/inference/graph_blockmodel_multicanonical.hh

@@ -40,6 +40,8 @@ using namespace std;
     ((S_min, , double, 0))                                                     \
     ((S_max, , double, 0))                                                     \
     ((f, , double, 0))                                                         \
+    ((time, , double, 0))                                                      \
+    ((refine, , bool, 0))                                                      \
     ((S, , double, 0))                                                         \
     ((E,, size_t, 0))                                                          \
     ((vlist,&, std::vector<size_t>&, 0))                                       \

src/graph/inference/graph_inference.cc

@@ -26,7 +26,7 @@ using namespace boost;
 using namespace graph_tool;
 
 void collect_vertex_marginals(GraphInterface& gi, boost::any ob,
-                              boost::any op)
+                              boost::any op, double update)
 {
     typedef vprop_map_t<int32_t>::type vmap_t;
     auto b = any_cast<vmap_t>(ob).get_unchecked();
@@ -34,22 +34,24 @@ void collect_vertex_marginals(GraphInterface& gi, boost::any ob,
     run_action<>()
         (gi, [&](auto& g, auto p)
          {
+             typename property_traits<decltype(p)>::value_type::value_type
+                 up = update;
              parallel_vertex_loop
                  (g,
                   [&](auto v)
                   {
-                     auto r = b[v];
-                     auto& pv = p[v];
-                     if (pv.size() <= size_t(r))
-                         pv.resize(r + 1);
-                     pv[r]++;
+                      auto r = b[v];
+                      auto& pv = p[v];
+                      if (pv.size() <= size_t(r))
+                          pv.resize(r + 1);
+                      pv[r] += up;
                   });
          },
          vertex_scalar_vector_properties())(op);
 }
 
 void collect_edge_marginals(GraphInterface& gi, size_t B, boost::any ob,
-                            boost::any op)
+                            boost::any op, double update)
 {
     typedef vprop_map_t<int32_t>::type vmap_t;
     auto b = any_cast<vmap_t>(ob).get_unchecked();
@@ -58,6 +60,8 @@ void collect_edge_marginals(GraphInterface& gi, size_t B, boost::any ob,
         (gi,
          [&](auto& g, auto p)
          {
+             typename property_traits<decltype(p)>::value_type::value_type
+                 up = update;
              parallel_edge_loop
                  (g,
                   [&](const auto& e)
@@ -72,7 +76,7 @@ void collect_edge_marginals(GraphInterface& gi, size_t B, boost::any ob,
                       if (pv.size() < B * B)
                           pv.resize(B * B);
                       size_t j = r + B * s;
-                      pv[j]++;
+                      pv[j] += up;
                   });
          },
          edge_scalar_vector_properties())(op);

src/graph/inference/multicanonical_loop.hh

@@ -49,7 +49,6 @@ auto multicanonical_sweep(MulticanonicalState& state, RNG& rng)
     int M = hist.size();
     double S_min = state._S_min;
     double S_max = state._S_max;
-    double f = state._f;
 
     auto get_bin = [&](double x) -> int
         {
@@ -100,7 +99,11 @@ auto multicanonical_sweep(MulticanonicalState& state, RNG& rng)
         }
 
         hist[i]++;
-        dens[i] += f;
+        dens[i] += state._f;
+
+        state._time += 1./M;
+        if (state._refine)
+            state._f = 1. / state._time;
     }
     return make_pair(S, nmoves);
 }

src/graph_tool/__init__.py

@@ -3086,7 +3086,7 @@ def _set_array_view(self, v):
     self.get_array()[:] = v
 
 vector_types = [Vector_bool, Vector_int16_t, Vector_int32_t, Vector_int64_t,
-                Vector_double, Vector_long_double]
+                Vector_double, Vector_long_double, Vector_size_t]
 for vt in vector_types:
     vt.a = property(_get_array_view, _set_array_view,
                     doc=r"""Shortcut to the `get_array` method as an attribute.""")

src/graph_tool/inference/blockmodel.py

@@ -286,7 +286,6 @@ class BlockState(object):
     def __setstate__(self, state):
         conv_pickle_state(state)
         self.__init__(**state)
-        return state
 
     def get_block_state(self, b=None, vweight=False, deg_corr=False, **kwargs):
         r"""Returns a :class:`~graph_tool.community.BlockState`` corresponding to the
@@ -950,9 +949,9 @@ class BlockState(object):
         return libinference.multicanonical_sweep(multicanonical_state,
                                                  self._state, _get_rng())
 
-    def multicanonical_sweep(self, m_state, f=1., c=1., niter=1,
-                             entropy_args={}, allow_empty=True, vertices=None,
-                             block_list=None, verbose=False):
+    def multicanonical_sweep(self, m_state, c=1., niter=1, entropy_args={},
+                             allow_empty=True, vertices=None, block_list=None,
+                             verbose=False):
         r"""Perform ``niter`` sweeps of a non-markovian multicanonical sampling using
         the Wang-Landau algorithm.
 
@@ -961,8 +960,6 @@ class BlockState(object):
         m_state : :class:`~graph_tool.inference.MulticanonicalState`
             :class:`~graph_tool.inference.MulticanonicalState` instance
             containing the current state of the Wang-Landau run.
-        f : ``float`` (optional, default: ``1.``)
-            Density of states update step.
         c : ``float`` (optional, default: ``1.``)
             Sampling parameter ``c`` for move proposals: For :math:`c\to 0` the
             blocks are sampled according to the local neighbourhood of a given
@@ -992,13 +989,18 @@ class BlockState(object):
             Entropy difference after the sweeps.
         nmoves : ``int``
             Number of vertices moved.
-
+ 
         References
         ----------
         .. [wang-efficient-2001] Fugao Wang, D. P. Landau, "An efficient, multiple
            range random walk algorithm to calculate the density of states", Phys.
            Rev. Lett. 86, 2050 (2001), :doi:`10.1103/PhysRevLett.86.2050`,
            :arxiv:`cond-mat/0011174`
+
+        .. [belardinelli-wang-2007] R. E. Belardinelli, V. D. Pereyra,
+           "Wang-Landau algorithm: A theoretical analysis of the saturation of
+           the error", J. Chem. Phys. 127, 184105 (2007),
+           :doi:`10.1063/1.2803061`, :arxiv:`cond-mat/0702414`
         """
 
         niter *= self.g.num_vertices()
@@ -1032,12 +1034,19 @@ class BlockState(object):
                                              ["xi_fast", "deg_dl_alt"]))
         multi_state.state = self._state
 
+        multi_state.f = m_state._f
+        multi_state.time = m_state._time
+        multi_state.refine = m_state._refine
         multi_state.S_min = m_state._S_min
         multi_state.S_max = m_state._S_max
         multi_state.hist = m_state._hist
         multi_state.dens = m_state._density
 
-        S, nmoves = self._multicanonical_sweep_dispatch(multi_state)
+        S, nmoves, f, time = \
+                self._multicanonical_sweep_dispatch(multi_state)
+
+        m_state._f = f
+        m_state._time = time
 
         if _bm_test():
             assert self._check_clabel(), "invalid clabel after sweep"
@@ -1154,7 +1163,7 @@ class BlockState(object):
             assert nB == B, "wrong number of groups after shrink: %d (should be %d)" % (nB, B)
         return state
 
-    def collect_edge_marginals(self, p=None):
+    def collect_edge_marginals(self, p=None, update=1.):
         r"""Collect the edge marginal histogram, which counts the number of times
         the endpoints of each node have been assigned to a given block pair.
 
@@ -1168,7 +1177,10 @@ class BlockState(object):
             membership counts.  Each vector entry corresponds to ``b[i] + B *
             b[j]``, where ``b`` is the block membership and ``i = min(source(e),
             target(e))`` and ``j = max(source(e), target(e))``. If not provided, an
-            empty histogram will be created.
+            empty histogram will be created
+        update : float (optional, default: ``1.``)
+            Each call increases the current count by the amount given by this
+            parameter.
 
         Returns
         -------
@@ -1199,15 +1211,16 @@ class BlockState(object):
         """
 
         if p is None:
-            p = self.g.new_ep("vector<int>")
+            p = self.g.new_ep("vector<double>")
 
         libinference.edge_marginals(self.g._Graph__graph,
                                     self.B,
                                     _prop("v", self.g, self.b),
-                                    _prop("e", self.g, p))
+                                    _prop("e", self.g, p),
+                                    update)
         return p
 
-    def collect_vertex_marginals(self, p=None):
+    def collect_vertex_marginals(self, p=None, update=1.):
         r"""Collect the vertex marginal histogram, which counts the number of times a
         node was assigned to a given block.
 
@@ -1219,6 +1232,9 @@ class BlockState(object):
         p : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
             Vertex property map with vector-type values, storing the previous block
             membership counts. If not provided, an empty histogram will be created.
+        update : float (optional, default: ``1.``)
+            Each call increases the current count by the amount given by this
+            parameter.
 
         Returns
         -------
@@ -1262,11 +1278,12 @@ class BlockState(object):
         """
 
         if p is None:
-            p = self.g.new_vp("vector<int>")
+            p = self.g.new_vp("vector<double>")
 
         libinference.vertex_marginals(self.g._Graph__graph,
                                       _prop("v", self.g, self.b),
-                                      _prop("v", self.g, p))
+                                      _prop("v", self.g, p),
+                                      update)
         return p
 
     def draw(self, **kwargs):

src/graph_tool/inference/blockmodel_em.py

@@ -122,6 +122,21 @@ class EMBlockState(object):
                     self.prs[r, s] = self.N * m[r, s] / (init_state.wr[r] * init_state.wr[s])
                     self.prs[s, r] = self.prs[r, s]
 
+    def __getstate__(self):
+        state = [self.g, self.B, self.vm, self.em_s, self.em_t, self.wr,
+                 self.prs]
+        return state
+
+    def __setstate__(self, state):
+        conv_pickle_state(state)
+        g, B, vm, em_s, em_t, wr, prs = state
+        self.__init__(g, B)
+        g.copy_property(vm, self.vm)
+        g.copy_property(em_s, self.em_s)
+        g.copy_property(em_t, self.em_t)
+        self.wr[:] = wr
+        self.prs[:,:] = prs
+
     def get_vertex_marginals(self):
         """Return the vertex marginals."""
         return self.vm

src/graph_tool/inference/layered_blockmodel.py

@@ -301,7 +301,6 @@ class LayeredBlockState(OverlapBlockState, BlockState):
     def __setstate__(self, state):
         conv_pickle_state(state)
         self.__init__(**state)
-        return state
 
     def __copy__(self):
         return self.copy()

src/graph_tool/inference/mcmc.py

@@ -381,6 +381,8 @@ class MulticanonicalState(object):
 
     Parameters
     ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be modelled.
     S_min : ``float``
         Minimum energy.
     S_max : ``float``
@@ -389,26 +391,34 @@ class MulticanonicalState(object):
         Number of bins.
     """
 
-    def __init__(self, S_min, S_max, nbins=1000):
+    def __init__(self, g, S_min, S_max, nbins=1000):
+        self._g = g
+        self._N = g.num_vertices()
         self._S_min = S_min
         self._S_max = S_max
         self._density = Vector_double()
         self._density.resize(nbins)
         self._hist = Vector_size_t()
         self._hist.resize(nbins)
+        self._perm_hist = numpy.zeros(nbins, dtype="int")
+        self._f = None
+        self._time = 0
+        self._refine = False
 
     def __getstate__(self):
-        state = [self._S_min, self._S_max,
-                 numpy.array(self._density.get_array()),
-                 numpy.array(self._hist.get_array())]
+        state = [self._g, self._S_min, self._S_max,
+                 numpy.array(self._density.a), numpy.array(self._hist.a),
+                 numpy.array(self._perm_hist), self._f, self._time,
+                 self._refine]
         return state
 
     def __setstate__(self, state):
-        S_min, S_max, density, hist = state
-        self.__init__(S_min, S_max, len(hist))
-        self._density.get_array()[:] = density
-        self._hist.get_array()[:] = hist
-        return state
+        g, S_min, S_max, density, hist, phist, self._f, self._time, \
+            self._refine = state
+        self.__init__(g, S_min, S_max, len(hist))
+        self._density.a[:] = density
+        self._hist.a[:] = hist
+        self._perm_hist[:] = phist
 
     def get_energies(self):
         "Get energy bounds."
@@ -416,7 +426,7 @@ class MulticanonicalState(object):
 
     def get_allowed_energies(self):
         "Get allowed energy bounds."
-        h = self._hist.get_array()
+        h = self._hist.a
         Ss = self.get_range()
         Ss = Ss[h > 0]
         return Ss[0], Ss[-1]
@@ -425,63 +435,66 @@ class MulticanonicalState(object):
         "Get energy range."
         return numpy.linspace(self._S_min, self._S_max, len(self._hist))
 
-    def get_density(self):
-        """Get density of states, normalized so that the **integral** over the energy
-        range is unity."""
-        r = numpy.array(self._density.get_array())
+    def get_density(self, B=None):
+        """Get density of states, normalized so that total sum is :math:`B^N`, where
+        :math:`B` is the number of groups, and :math:`N` is the number of
+        nodes. If not supplied :math:`B=N` is assumed.
+        """
+        r = numpy.array(self._density.a)
         r -= r.max()
-        N = len(r)
-        dS = (self._S_max - self._S_min) / N
-        I = exp(r).sum() * dS
-        r -= log(I)
+        r -= log(exp(r).sum())
+        if B == None:
+            B = self._g.num_vertices()
+        r += self._g.num_vertices() * log(B)
         return r
 
-    def get_prob(self, S):
-        r = self.get_density()
-        dS = (self._S_max - self._S_min) / N
+    def get_entropy(self, S, B=None):
+        r = self.get_density(B)
+        dS = (self._S_max - self._S_min) / len(r)
         j = round((S - self._S_min) / dS)
         return r[j]
 
     def get_hist(self):
         "Get energy histogram."
-        return numpy.array(self._hist.get_array())
+        return numpy.array(self._hist.a)
+
+    def get_perm_hist(self):
+        "Get permanent energy histogram."
+        return self._perm_hist
 
-    def get_flatness(self):
+    def get_flatness(self, use_ent=True, h=None):
         "Get energy histogram flatness."
-        h = self._hist.get_array()
+        if h is None:
+            h = self._hist.a
         if h.sum() == 0:
-            return numpy.inf
+            return 0
         h = array(h[h>0], dtype="float")
-        h /= h.sum()
-        S = -(h * log(h)).sum()
-        return len(h) / exp(S) - 1
-
-    def get_mean(self):
-        "Get energy mean."
-        r = self.get_density()
-        N = len(r)
-        dS = (self._S_max - self._S_min) / N
-        Ss = numpy.linspace(self._S_min, self._S_max, N)
-        return (Ss * exp(r) * dS).sum()
-
-    def get_posterior(self):
-        "Get posterior mean."
-        r = self.get_density()
-        N = len(r)
-        dS = (self._S_max - self._S_min) / N
-        Ss = numpy.linspace(self._S_min, self._S_max, N)
+        if not use_ent:
+            h_mean = h.mean()
+            return h.min() / h_mean
+        else:
+            h /= h.sum()
+            S = -(h * log(h)).sum()
+            return exp(S - log(len(h)))
+
+    def get_posterior(self, N=None):
+        "Get posterior probability."
+        r = self.get_density(N)
+        Ss = numpy.linspace(self._S_min, self._S_max, len(r))
         y = -Ss + r
         y_max = y.max()
         y -= y_max
-        return y_max + log((exp(y) * dS).sum())
+        return y_max + log(exp(y).sum())
 
     def reset_hist(self):
         "Reset energy histogram."
-        self._hist.get_array()[:] = 0
+        self._perm_hist += asarray(self._hist.a[:], dtype="int")
+        self._hist.a[:] = 0
 
-def multicanonical_equilibrate(state, m_state, f_range=(1., 1e-6), r=2,
-                               flatness=.01, callback=None,
-                               multicanonical_args={}, verbose=False):
+def multicanonical_equilibrate(state, m_state, f_range=(1., 1e-6),
+                               f_refine=True, r=2, flatness=.99, use_ent=True,
+                               callback=None, multicanonical_args={},
+                               verbose=False):
     r"""Equilibrate a multicanonical Monte Carlo sampling using the Wang-Landau
      algorithm.
 
@@ -493,11 +506,17 @@ def multicanonical_equilibrate(state, m_state, f_range=(1., 1e-6), r=2,
         Initial multicanonical state, where the state density will be stored.
     f_range : ``tuple`` of two floats (optional, default: ``(1., 1e-6)``)
         Range of density updates.
+    f_refine : ``bool`` (optional, default: ``True``)
+        If ``True``, the refinement steps described in [belardinelli-wang-2007]_
+        will be used.
     r : ``float`` (optional, default: ``2.``)
         Greediness of convergence. At each iteration, the density updates will
         be reduced by a factor ``r``.
-    flatness : ``float`` (optional, default: ``1e-3``)
+    flatness : ``float`` (optional, default: ``.99``)
         Sufficient histogram flatness threshold used to continue the algorithm.
+    use_ent : ``bool`` (optional, default: ``True``)
+        If ``True``, the histogram entropy will be used to determine flatness,
+        otherwise the smallest count relative to the mean will be used.
     callback : ``function`` (optional, default: ``None``)
         If given, this function will be called after each iteration. The
         function must accept the current ``state`` and ``m_state`` as arguments.
@@ -524,24 +543,35 @@ def multicanonical_equilibrate(state, m_state, f_range=(1., 1e-6), r=2,
        range random walk algorithm to calculate the density of states", Phys.
        Rev. Lett. 86, 2050 (2001), :doi:`10.1103/PhysRevLett.86.2050`,
        :arxiv:`cond-mat/0011174`
+    .. [belardinelli-wang-2007] R. E. Belardinelli, V. D. Pereyra,
+       "Wang-Landau algorithm: A theoretical analysis of the saturation of
+       the error", J. Chem. Phys. 127, 184105 (2007),
+       :doi:`10.1063/1.2803061`, :arxiv:`cond-mat/0702414`
     """
 
     count = 0
-    f = f_range[0]
-    while f >= f_range[1]:
-        state.multicanonical_sweep(m_state, **overlay(multicanonical_args, f=f))
-        hf = m_state.get_flatness()
-        if hf < flatness:
-            f /= r
-            if f >= f_range[1]:
-                m_state.reset_hist()
+    if m_state._f is None:
+        m_state._f = f_range[0]
+    while m_state._f >= f_range[1]:
+        state.multicanonical_sweep(m_state, **multicanonical_args)
+        hf = m_state.get_flatness(use_ent=use_ent)
 
         if callback is not None:
-            calback(state, m_state)
+            callback(state, m_state)
 
-        count += 1
         if check_verbose(verbose):
             print(verbose_pad(verbose) +
-                  "iter: %d  f: %g  flatness: %g" % (count, f, hf))
+                  "count: %d  time: %#8.8g  f: %#8.8g  flatness: %#8.8g  S: %#8.8g" % \
+                  (count, m_state._time, m_state._f, hf,
+                   state.entropy(multicanonical_args.get("entropy_args", {}))))
+
+        if not m_state._refine:
+            if hf > flatness:
+                m_state._f /= r
+                if m_state._f >= f_range[1]:
+                    m_state.reset_hist()
+                if f_refine and m_state._f < 1 / m_state._time:
+                    m_state._refine = True
+        count += 1
 
     return count
\ No newline at end of file

src/graph_tool/inference/nested_blockmodel.py

@@ -101,7 +101,6 @@ class NestedBlockState(object):
     def __setstate__(self, state):
         conv_pickle_state(state)
         self.__init__(**overlay(dmask(state, ["kwargs"]), **state["kwargs"]))
-        return state
 
     def project_partition(self, j, l):
         """Project partition of level ``j`` onto level ``l``, and return it."""

src/graph_tool/inference/overlap_blockmodel.py

@@ -259,7 +259,6 @@ class OverlapBlockState(BlockState):
     def __setstate__(self, state):
         conv_pickle_state(state)
         self.__init__(**state)
-        return state
 
     def get_block_state(self, b=None, vweight=False, overlap=False,
                         deg_corr=False, **kwargs):