From bc4b3707a15385849e3362763f408049a4cb56cf Mon Sep 17 00:00:00 2001 From: Tiago de Paula Peixoto Date: Tue, 21 May 2013 19:47:40 +0200 Subject: [PATCH] Implement fully random moves in blockmodel.py --- src/graph/community/graph_blockmodel.cc | 21 +++--- src/graph/community/graph_blockmodel.hh | 10 +-- src/graph_tool/community/blockmodel.py | 88 +++++++++++++++---------- 3 files changed, 69 insertions(+), 50 deletions(-) diff --git a/src/graph/community/graph_blockmodel.cc b/src/graph/community/graph_blockmodel.cc index 54654b2c..e1e8374f 100644 --- a/src/graph/community/graph_blockmodel.cc +++ b/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& 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::type @@ -484,8 +487,8 @@ python::object do_move_sweep(GraphInterface& gi, GraphInterface& bgi, (gi, boost::bind(move_sweep_dispatch (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); } diff --git a/src/graph/community/graph_blockmodel.hh b/src/graph/community/graph_blockmodel.hh index ba3a9e2f..3815fb50 100644 --- a/src/graph/community/graph_blockmodel.hh +++ b/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& 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::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 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); diff --git a/src/graph_tool/community/blockmodel.py b/src/graph_tool/community/blockmodel.py index 8ac695e0..7e8d7b29 100644 --- a/src/graph_tool/community/blockmodel.py +++ b/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 -- GitLab