Implement fully random moves in blockmodel.py

src/graph/community/graph_blockmodel.cc

@@ -370,13 +370,14 @@ struct move_sweep_dispatch
     move_sweep_dispatch(Eprop eweight, Vprop vweight, boost::any egroups,
                         VEprop esrcpos, VEprop etgtpos, Vprop label, size_t L,
                         vector<int>& vlist, bool deg_corr, double beta,
-                        bool sequential, bool verbose, size_t max_edge_index,
-                        rng_t& rng, double& S, size_t& nmoves,
-                        GraphInterface& bgi)
+                        bool sequential, bool random_move, bool verbose,
+                        size_t max_edge_index, rng_t& rng, double& S,
+                        size_t& nmoves, GraphInterface& bgi)
 
         : eweight(eweight), vweight(vweight), oegroups(egroups), esrcpos(esrcpos),
           etgtpos(etgtpos), label(label), L(L), vlist(vlist),
-          deg_corr(deg_corr), beta(beta), sequential(sequential), verbose(verbose),
+          deg_corr(deg_corr), beta(beta), sequential(sequential),
+          random_move(random_move), verbose(verbose),
           max_edge_index(max_edge_index), rng(rng), S(S), nmoves(nmoves), bgi(bgi)
     {}
 
@@ -392,6 +393,7 @@ struct move_sweep_dispatch
     bool deg_corr;
     double beta;
     bool sequential;
+    bool random_move;
     bool verbose;
     size_t max_edge_index;
     rng_t& rng;
@@ -441,7 +443,8 @@ struct move_sweep_dispatch
                    egroups.get_unchecked(num_vertices(bg)),
                    esrcpos.get_unchecked(max_edge_index + 1),
                    etgtpos.get_unchecked(max_edge_index + 1),
-                   g, bg, emat, sequential, verbose, rng, S, nmoves);
+                   g, bg, emat, sequential, random_move, verbose,
+                   rng, S, nmoves);
     }
 
 };
@@ -454,8 +457,8 @@ python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi,
                              bool deg_corr, boost::any oeweight,
                              boost::any ovweight, boost::any oegroups,
                              boost::any oesrcpos, boost::any oetgtpos,
-                             double beta, bool sequential, bool verbose,
-                             rng_t& rng)
+                             double beta, bool sequential, bool random_move,
+                             bool verbose, rng_t& rng)
 {
     typedef property_map_type::apply<int32_t,
                                      GraphInterface::vertex_index_map_t>::type
@@ -484,8 +487,8 @@ python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi,
         (gi, boost::bind<void>(move_sweep_dispatch<emap_t, vmap_t, vemap_t>
                                (eweight, vweight, oegroups, esrcpos, etgtpos,
                                 label, L, vlist, deg_corr, beta,
-                                sequential, verbose, gi.GetMaxEdgeIndex(),
-                                rng, S, nmoves, bgi),
+                                sequential, random_move,  verbose,
+                                gi.GetMaxEdgeIndex(), rng, S, nmoves, bgi),
                                mrs, mrp, mrm, wr, b, _1, ref(emat)))();
     return python::make_tuple(S, nmoves);
 }

src/graph/community/graph_blockmodel.hh

@@ -799,8 +799,8 @@ void move_sweep(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Vprop label,
                 size_t L, vector<int>& vlist, bool deg_corr, double beta,
                 Eprop eweight, Vprop vweight, EVprop egroups, VEprop esrcpos,
                 VEprop etgtpos, Graph& g, BGraph& bg, EMat& emat,
-                bool sequential, bool verbose, RNG& rng, double& S,
-                size_t& nmoves)
+                bool sequential, bool random_move, bool verbose, RNG& rng,
+                double& S, size_t& nmoves)
 {
     typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
     nmoves = 0;
@@ -841,7 +841,7 @@ void move_sweep(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Vprop label,
         else
             p_rand = B / double(mrp[t] + B);
 
-        if (rand_real() >= p_rand)
+        if (rand_real() >= p_rand && !random_move)
         {
             std::tr1::uniform_int<size_t> urand(0, egroups[t].size() - 1);
 
@@ -857,12 +857,12 @@ void move_sweep(Eprop mrs, Vprop mrp, Vprop mrm, Vprop wr, Vprop b, Vprop label,
         if (s == r)
             continue;
 
-        double pf = isinf(beta) ? 1 : get_move_prob(v, s, b, mrs, mrp, mrm, emat, g, B);
+        double pf = (isinf(beta) || random_move) ? 1 : get_move_prob(v, s, b, mrs, mrp, mrm, emat, g, B);
 
         double dS = move_vertex(v, s, mrs, mrp, mrm, wr, b, deg_corr, eweight,
                                 vweight, g, bg, emat, true);
 
-        double pb = isinf(beta) ? 1 : get_move_prob(v, r, b, mrs, mrp, mrm, emat, g, B);
+        double pb = (isinf(beta) || random_move) ? 1 : get_move_prob(v, r, b, mrs, mrp, mrm, emat, g, B);
 
         double a = isinf(beta) ? -dS : -beta * dS + log(pb) - log(pf);
 

src/graph_tool/community/blockmodel.py

@@ -716,7 +716,7 @@ def min_dist(state, n=0):
     return min_d, r, s
 
 
-def mcmc_sweep(state, beta=1., sequential=True, verbose=False, vertices=None):
+def mcmc_sweep(state, beta=1., sequential=True, vertices=None, random_move=False, verbose=False):
     r"""Performs a Monte Carlo Markov chain sweep on the network, to sample the block partition according to a probability :math:`\propto e^{-\beta \mathcal{S}_{t/c}}`, where :math:`\mathcal{S}_{t/c}` is the blockmodel entropy.
 
     Parameters
@@ -730,6 +730,14 @@ def mcmc_sweep(state, beta=1., sequential=True, verbose=False, vertices=None):
         random order. Otherwise a total of `N` moves attempts are made, where
         `N` is the number of vertices, where each vertex can be selected with
         equal probability.
+    random_move : ``bool`` (optional, default: ``False``)
+        If ``True``, the proposed moves will attempt to place the vertices in
+        fully randomly-chosen blocks. If ``False``, the proposed moves will be
+        chosen with a probability depending on the membership of the neighbours
+        and the currently-inferred block structure.
+    vertices: ``list of ints`` (optional, default: ``None``)
+        A list of vertices which will be attempted to be moved. If ``None``, all
+        vertices will be attempted.
     verbose : ``bool`` (optional, default: ``False``)
         If ``True``, verbose information is displayed.
 
@@ -838,6 +846,10 @@ def mcmc_sweep(state, beta=1., sequential=True, verbose=False, vertices=None):
     if vertices is None:
         vertices = libcommunity.get_vector(state.g.num_vertices())
         vertices.a = state.g.vertex_index.copy("int").fa
+    else:
+        vlist = libcommunity.get_vector(len(vertices))
+        vlist.a = vertices
+        vertices = vlist
 
     dS, nmoves = libcommunity.move_sweep(state.g._Graph__graph,
                                          state.bg._Graph__graph,
@@ -856,62 +868,66 @@ def mcmc_sweep(state, beta=1., sequential=True, verbose=False, vertices=None):
                                          state.egroups,
                                          _prop("e", state.g, state.esrcpos),
                                          _prop("e", state.g, state.etgtpos),
-                                         float(beta), sequential,
+                                         float(beta), sequential, random_move,
                                          verbose, _get_rng())
     return dS / state.E, nmoves
 
 
 
 
-def mc_get_dl(state, nsweep, greedy, rng, checkpoint=None, anneal=1,
-              verbose=False):
+def mc_get_dl(state, nsweep, greedy, rng, checkpoint, checkpoint_state,
+              anneal=1, verbose=False):
 
     if len(state.vertices) == 1:
         return state._BlockState__min_dl()
 
-    S = state._BlockState__min_dl()
+    B = len(state.vertices)
+    if checkpoint_state is None:
+        checkpoint_state = defaultdict(dict)
+    if "S" in checkpoint_state[B]:
+        S = checkpoint_state[B]["S"]
+        niter = checkpoint_state[B]["niter"]
+    else:
+        S = state._BlockState__min_dl()
+        checkpoint_state[B]["S"] = S
+        niter = 0
+        checkpoint_state[B]["niter"] = niter
 
     if nsweep >= 0:
-        for i in range(nsweep):
-            delta, nmoves = mcmc_sweep(state, beta=1, rng=rng)
-            if checkpoint is not None:
-                checkpoint(state, S, delta, nmoves)
+        ntotal = nsweep if not greedy else 2 * nsweep
+        for i in range(ntotal):
+            if i < niter:
+                continue
+            if i < nsweep:
+                beta = 1
+            else:
+                beta = float("inf") # greedy
+            delta, nmoves = mcmc_sweep(state, beta=beta, rng=rng)
             S += delta
-        if greedy:
-            for i in range(nsweep):
-                delta, nmoves = mcmc_sweep(state, beta=float("inf"))
-                if checkpoint is not None:
-                    checkpoint(state, S, delta, nmoves)
-                    S += delta
+            niter += 1
+            checkpoint_state[B]["S"] = S
+            checkpoint_state[B]["niter"] = niter
+            if checkpoint is not None:
+                checkpoint(state, S, delta, nmoves, checkpoint_state)
     else:
         # adaptive mode
-        min_dl = S
-        max_dl = S
-        count = 0
-        time = 0
-        bump = False
+        min_dl = checkpoint_state[B].get("min_dl", S)
+        max_dl = checkpoint_state[B].get("max_dl", S)
+        count = checkpoint_state[B].get("count", 0)
+        time = checkpoint_state[B].get("time", 0)
+        bump = checkpoint_state[B].get("bump", False)
 
-        beta = 1.
+        beta =  checkpoint_state[B].get("beta", 1.)
 
-        last_min = min_dl
+        last_min = checkpoint_state[B].get("last_min", min_dl)
+
+        greedy = checkpoint_state[B].get("greedy", False)
 
         if verbose:
             print("beta = %g" % beta)
         while True:
-            delta, nmoves = mcmc_sweep(state, beta=beta)
-            if checkpoint is not None:
-                checkpoint(state, S, delta, nmoves)
-            S += delta
-
-            if S < min_dl:
-                min_dl = S
-                count = 0
-            elif S > max_dl:
-                max_dl = S
-                count = 0
-            else:
-                count += 1
-
+            if greedy:
+                break
             if count > abs(nsweep):
                 if not bump:
                     min_dl = max_dl = S