blockmodel.py 107 KB
Newer Older
1
2
3
4
5
#! /usr/bin/env python
# -*- coding: utf-8 -*-
#
# graph_tool -- a general graph manipulation python module
#
Tiago Peixoto's avatar
Tiago Peixoto committed
6
# Copyright (C) 2006-2014 Tiago de Paula Peixoto <tiago@skewed.de>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#
# 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

26
27
from .. import _degree, _prop, Graph, GraphView, libcore, _get_rng, PropertyMap
from .. stats import label_self_loops
28
from .. spectral import adjacency
29
30
import random
from numpy import *
31
import numpy
32
33
from scipy.optimize import fsolve, fminbound
import scipy.special
34
from collections import defaultdict
35
36
import copy
import heapq
37
38
39
40

from .. dl_import import dl_import
dl_import("from . import libgraph_tool_community as libcommunity")

41
__test__ = False
42

43
44
45
46
47
48
49
50
51
52
53
54
def get_block_graph(g, B, b, vcount, ecount):
    cg, br, vcount, ecount = condensation_graph(g, b,
                                                vweight=vcount,
                                                eweight=ecount,
                                                self_loops=True)[:4]
    cg.vp["count"] = vcount
    cg.ep["count"] = ecount
    cg = Graph(cg, vorder=br)

    cg.add_vertex(B - cg.num_vertices())
    return cg

55
56
57
58
59
60
61
62
class BlockState(object):
    r"""This class encapsulates the block state of a given graph.

    This must be instantiated and used by functions such as :func:`mcmc_sweep`.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
63
        Graph to be modelled.
64
    eweight : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
65
        Edge multiplicities (for multigraphs or block graphs).
66
    vweight : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
67
        Vertex multiplicities (for block graphs).
68
69
70
71
72
73
74
75
    b : :class:`~graph_tool.PropertyMap` (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. If not supplied it will be either obtained from the
        parameter ``b``, or set to the maximum possible value according to the
        minimum description length.
    clabel : :class:`~graph_tool.PropertyMap` (optional, default: ``None``)
76
77
        Constraint labels on the vertices. If supplied, vertices with different
        label values will not be clustered in the same group.
78
79
80
    deg_corr : ``bool`` (optional, default: ``True``)
        If ``True``, the degree-corrected version of the blockmodel ensemble will
        be assumed, otherwise the traditional variant will be used.
81
82
83
84
    max_BE : ``int`` (optional, default: ``1000``)
        If the number of blocks exceeds this number, a sparse representation of
        the block graph is used, which is slightly less efficient, but uses less
        memory,
85
86
    """

87
88
    _state_ref_count = 0

89
    def __init__(self, g, eweight=None, vweight=None, b=None,
90
91
92
                 B=None, clabel=None, deg_corr=True,
                 max_BE=1000, **kwargs):

93
94
        BlockState._state_ref_count += 1

95
        # initialize weights to unity, if necessary
96
97
        if eweight is None:
            eweight = g.new_edge_property("int")
98
            eweight.fa = 1
99
100
101
102
        elif eweight.value_type() != "int32_t":
            eweight = eweight.copy(value_type="int32_t")
        if vweight is None:
            vweight = g.new_vertex_property("int")
103
            vweight.fa = 1
104
105
        elif vweight.value_type() != "int32_t":
            vweight = vweight.copy(value_type="int32_t")
106
107
108
109
110
111
112
        self.eweight = g.own_property(eweight)
        self.vweight = g.own_property(vweight)

        self.is_weighted = True if self.eweight.fa.max() > 1 else False

        # configure the main graph and block model parameters
        self.g = g
113

114
115
        self.E = int(self.eweight.fa.sum())
        self.N = int(self.vweight.fa.sum())
116
117
118

        self.deg_corr = deg_corr

119
120
121
122
        # ensure we have at most as many blocks as nodes
        if B is not None:
            B = min(B, self.g.num_vertices())

123
        if b is None:
124
            # create a random partition into B blocks.
125
126
            if B is None:
                B = get_max_B(self.N, self.E, directed=g.is_directed())
127
128
            B = min(B, self.g.num_vertices())
            ba = random.randint(0, B, self.g.num_vertices())
129
            ba[:B] = arange(B)        # avoid empty blocks
130
131
            if B < self.g.num_vertices():
                random.shuffle(ba)
132
            b = g.new_vertex_property("int")
133
            b.fa = ba
134
135
            self.b = b
        else:
136
137
138
139
140
141
            # if a partition is available, we will incorporate it.
            if isinstance(b, numpy.ndarray):
                self.b = g.new_vertex_property("int")
                self.b.fa = b
            else:
                self.b = b = g.own_property(b.copy(value_type="int"))
142
            if B is None:
143
144
145
146
                B = int(self.b.fa.max()) + 1

        # if B > self.N:
        #     raise ValueError("B > N!")
147

148
        if self.b.fa.max() >= B:
149
150
151
            raise ValueError("Maximum value of b is larger or equal to B!")

        # Construct block-graph
152
        self.bg = get_block_graph(g, B, self.b, self.vweight, self.eweight)
153
154
155
156
        self.bg.set_fast_edge_removal()

        self.mrs = self.bg.ep["count"]
        self.wr = self.bg.vp["count"]
157

158
159
160
161
162
163
164
165
166
        del self.bg.ep["count"]
        del self.bg.vp["count"]

        self.mrp = self.bg.degree_property_map("out", weight=self.mrs)

        if g.is_directed():
            self.mrm = self.bg.degree_property_map("in", weight=self.mrs)
        else:
            self.mrm = self.mrp
167
168
169

        self.vertices = libcommunity.get_vector(B)
        self.vertices.a = arange(B)
170
        self.B = B
171

172
173
174
175
176
177
178
        if clabel is not None:
            if isinstance(clabel, PropertyMap):
                self.clabel = clabel
            else:
                self.clabel = self.g.new_vertex_property("int")
                self.clabel.a = clabel
        else:
179
180
181
182
183
184
185
            self.clabel = self.g.new_vertex_property("int")

        self.emat = None
        if max_BE is None:
            max_BE = 1000
        self.max_BE = max_BE

186
187
        self.overlap = False

188
189
190
191
        # used by mcmc_sweep()
        self.egroups = None
        self.nsampler = None
        self.sweep_vertices = None
192
193
        self.overlap_stats = libcommunity.overlap_stats()
        self.partition_stats = libcommunity.partition_stats()
194

195
196
197
        # computation cache
        libcommunity.init_safelog(int(5 * max(self.E, self.N)))
        libcommunity.init_xlogx(int(5 * max(self.E, self.N)))
198
        libcommunity.init_lgamma(int(3 * max(self.E, self.N)))
199

200
    def __del__(self):
Tiago Peixoto's avatar
Tiago Peixoto committed
201
202
203
204
205
206
207
208
        try:
            BlockState._state_ref_count -= 1
            if BlockState._state_ref_count == 0:
                libcommunity.clear_safelog()
                libcommunity.clear_xlogx()
                libcommunity.clear_lgamma()
        except (ValueError, AttributeError):
            pass
209

210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
    def __repr__(self):
        return "<BlockState object with %d blocks,%s for graph %s, at 0x%x>" % \
            (self.B, " degree corrected," if self.deg_corr else "", str(self.g),
             id(self))


    def __init_partition_stats(self, empty=True):
        if not empty:
            self.partition_stats = libcommunity.init_partition_stats(self.g._Graph__graph,
                                                                     _prop("v", self.g, self.b),
                                                                     _prop("e", self.g, self.eweight),
                                                                     self.N, self.B)
        else:
            self.partition_stats = libcommunity.partition_stats()



    def copy(self, b=None, B=None, deg_corr=None, clabel=None, overlap=False):
        r"""Copies the block state. The parameters override the state properties, and
         have the same meaning as in the constructor. If ``overlap=True`` an
         instance of :class:`~graph_tool.community.OverlapBlockState` is
         returned."""

        if not overlap:
            state = BlockState(self.g,
                               eweight=self.eweight,
                               vweight=self.vweight,
                               b=self.b.copy() if b is None else b,
                               B=(self.B if b is None else None) if B is None else B,
                               clabel=self.clabel if clabel is None else clabel,
                               deg_corr=self.deg_corr if deg_corr is None else deg_corr,
                               max_BE=self.max_BE)
        else:
            state = OverlapBlockState(self.g,
                                      b=b if b is not None else self.b,
                                      B=(self.B if b is None else None) if B is None else B,
                                      clabel=self.clabel if clabel is None else clabel,
                                      deg_corr=self.deg_corr if deg_corr is None else deg_corr,
                                      max_BE=self.max_BE)

        if not state.__check_clabel():
            b = state.b.a + state.clabel.a * state.B
            continuous_map(b)
            state = state.copy(b=b)

            if __test__:
                assert state.__check_clabel()

        return state


    def __getstate__(self):
        state = dict(g=self.g,
                     eweight=self.eweight,
                     vweight=self.vweight,
                     b=self.b,
                     B=self.B,
                     clabel=self.clabel,
                     deg_corr=self.deg_corr,
                     max_BE=self.max_BE)
        return state

    def __setstate__(self, state):
        self.__init__(**state)
        return state

    def get_block_state(self, b=None, vweight=False, deg_corr=False, overlap=False):
        r"""Returns a :class:`~graph_tool.community.BlockState`` corresponding to the
        block graph. The parameters have the same meaning as the in the constructor."""


        state = BlockState(self.bg, eweight=self.mrs,
                           vweight=self.wr if vweight else None,
                           b=self.bg.vertex_index.copy("int") if b is None else b,
                           clabel=self.get_bclabel(),
                           deg_corr=deg_corr,
                           max_BE=self.max_BE)
        if overlap:
            state = state.copy(overlap=True)
        n_map = self.b.copy()
        return state, n_map

    def get_bclabel(self):
        r"""Returns a :class:`~graph_tool.PropertyMap`` corresponding to constraint
        labels for the block graph."""

        bclabel = self.bg.new_vertex_property("int")
        reverse_map(self.b, bclabel)
        pmap(bclabel, self.clabel)
        return bclabel

    def __check_clabel(self):
        b = self.b.a + self.clabel.a * self.B
        continuous_map(b)
        b2 = self.b.copy()
        continuous_map(b2.a)
        return (b == b2.a).all()

308
309
310
311
    def __get_emat(self):
        if self.emat is None:
            self.__regen_emat()
        return self.emat
312
313

    def __regen_emat(self):
314
315
316
317
        if self.B <= self.max_BE:
            self.emat = libcommunity.create_emat(self.bg._Graph__graph)
        else:
            self.emat = libcommunity.create_ehash(self.bg._Graph__graph)
318

319
    def __build_egroups(self, empty=False):
320
321
        self.esrcpos = self.g.new_edge_property("int")
        self.etgtpos = self.g.new_edge_property("int")
322

323
        self.egroups = libcommunity.build_egroups(self.g._Graph__graph,
324
325
326
327
328
                                                  self.bg._Graph__graph,
                                                  _prop("v", self.g, self.b),
                                                  _prop("e", self.g, self.eweight),
                                                  _prop("e", self.g, self.esrcpos),
                                                  _prop("e", self.g, self.etgtpos),
329
                                                  self.is_weighted, empty)
330

331
    def __build_nsampler(self, empty=False):
332
        self.nsampler = libcommunity.init_neighbour_sampler(self.g._Graph__graph,
333
                                                            _prop("e", self.g, self.eweight),
334
                                                            True, empty)
335
336
337
338
339
340

    def __cleanup_bg(self):
        emask = self.bg.new_edge_property("bool")
        emask.a = self.mrs.a[:len(emask.a)] > 0
        self.bg.set_edge_filter(emask)
        self.bg.purge_edges()
341
        self.emat = None
342
343
344
345
346
347
348
349
350
351

    def get_blocks(self):
        r"""Returns the property map which contains the block labels for each vertex."""
        return self.b

    def get_bg(self):
        r"""Returns the block graph."""
        return self.bg

    def get_ers(self):
352
353
        r"""Returns the edge property map of the block graph which contains the :math:`e_{rs}` matrix entries.
        For undirected graphs, the diagonal values (self-loops) contain :math:`e_{rr}/2`."""
354
355
356
357
358
359
360
361
        return self.mrs

    def get_er(self):
        r"""Returns the vertex property map of the block graph which contains the number
        :math:`e_r` of half-edges incident on block :math:`r`. If the graph is
        directed, a pair of property maps is returned, with the number of
        out-edges :math:`e^+_r` and in-edges :math:`e^-_r`, respectively."""
        if self.bg.is_directed():
362
            return self.mrp, self.mrm
363
364
365
366
367
368
369
        else:
            return self.mrp

    def get_nr(self):
        r"""Returns the vertex property map of the block graph which contains the block sizes :math:`n_r`."""
        return self.wr

370
371
372
    def entropy(self, complete=True, dl=False, partition_dl=False, dense=False,
                multigraph=True, norm=True, dl_ent=False, **kwargs):
        r"""Calculate the entropy associated with the current block partition.
373
374
375
376
377
378
379
380

        Parameters
        ----------
        complete : ``bool`` (optional, default: ``False``)
            If ``True``, the complete entropy will be returned, including constant
            terms not relevant to the block partition.
        dl : ``bool`` (optional, default: ``False``)
            If ``True``, the full description length will be returned.
381
382
383
        partition_dl : ``bool`` (optional, default: ``False``)
            If ``True``, and ``dl == True`` only the partition description
            length will be considered.
384
385
386
        dense : ``bool`` (optional, default: ``False``)
            If ``True``, the "dense" variant of the entropy will be computed.
        multigraph : ``bool`` (optional, default: ``False``)
387
388
389
390
391
392
393
            If ``True``, the multigraph entropy will be used.
        norm : ``bool`` (optional, default: ``True``)
            If ``True``, the entropy will be "normalized" by dividing by the
            number of edges.
        dl_ent : ``bool`` (optional, default: ``False``)
            If ``True``, the description length of the degree sequence will be
            approximated by its entropy.
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419

        Notes
        -----
        For the traditional blockmodel (``deg_corr == False``), the entropy is
        given by

        .. math::

          \mathcal{S}_t &\cong E - \frac{1}{2} \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{n_rn_s}\right), \\
          \mathcal{S}^d_t &\cong E - \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{n_rn_s}\right),

        for undirected and directed graphs, respectively, where :math:`e_{rs}`
        is the number of edges from block :math:`r` to :math:`s` (or the number
        of half-edges for the undirected case when :math:`r=s`), and :math:`n_r`
        is the number of vertices in block :math:`r` .

        For the degree-corrected variant with "hard" degree constraints the
        equivalent expressions are

        .. math::

            \mathcal{S}_c &\cong -E -\sum_kN_k\ln k! - \frac{1}{2} \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{e_re_s}\right), \\
            \mathcal{S}^d_c &\cong -E -\sum_{k^+}N_{k^+}\ln k^+!  -\sum_{k^-}N_{k^-}\ln k^-! - \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{e^+_re^-_s}\right),

        where :math:`e_r = \sum_se_{rs}` is the number of half-edges incident on
        block :math:`r`, and :math:`e^+_r = \sum_se_{rs}` and :math:`e^-_r =
420
        \sum_se_{sr}` are the numbers of out- and in-edges adjacent to block
421
422
        :math:`r`, respectively.

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
        If ``dense == False`` and ``multigraph == True``, the entropy used will
        be of the "Poisson" model, with the additional term:

        .. math::

            {\mathcal{S}_{cm}^{(d)}} = \mathcal{S}_c^{(d)} + \sum_{i>j} \ln A_{ij}! + \sum_i \ln A_{ii}!!


        If ``dl == True``, the description length :math:`\mathcal{L}_t` of the
        model will be returned as well, as described in
        :func:`model_entropy`. Note that for the degree-corrected version the
        description length is

        .. math::

            \mathcal{L}_c = \mathcal{L}_t + \sum_r\min\left(\mathcal{L}^{(1)}_r, \mathcal{L}^{(2)}_r\right),

        with

        .. math::

              \mathcal{L}^{(1)}_r &= \ln{\left(\!\!{n_r \choose e_r}\!\!\right)}, \\
              \mathcal{L}^{(2)}_r &= \ln\Xi_r + \ln n_r! - \sum_k \ln n^r_k!,

        and :math:`\ln\Xi_r \simeq 2\sqrt{\zeta(2)e_r}`, where :math:`\zeta(x)`
        is the `Riemann zeta function
        <https://en.wikipedia.org/wiki/Riemann_zeta_function>`_, and
        :math:`n^r_k` is the number of nodes in block :math:`r` with degree
        :math:`k`. For directed graphs we have instead :math:`k \to (k^-, k^+)`,
        and :math:`\ln\Xi_r \to \ln\Xi^+_r + \ln\Xi^-_r` with :math:`\ln\Xi_r^+
        \simeq 2\sqrt{\zeta(2)e^+_r}` and :math:`\ln\Xi_r^- \simeq
        2\sqrt{\zeta(2)e^-_r}`.

        If ``dl_ent=True`` is passed, this will be approximated instead by
457
458
459

        .. math::

460
            \mathcal{L}_c \simeq \mathcal{L}_t - \sum_rn_r\sum_kp^r_k\ln p^r_k,
461

462
        where :math:`p^r_k = n^r_k / n_r`.
463

464
465
        If the "dense" entropies are requested (``dense=True``), they will be
        computed as
466
467
468
469
470
471
472
473
474
475
476
477
478

        .. math::

            \mathcal{S}_t  &= \sum_{r>s} \ln{\textstyle {n_rn_s \choose e_{rs}}} + \sum_r \ln{\textstyle {{n_r\choose 2}\choose e_{rr}/2}}\\
            \mathcal{S}^d_t  &= \sum_{rs} \ln{\textstyle {n_rn_s \choose e_{rs}}},

        for simple graphs, and

        .. math::

            \mathcal{S}_m  &= \sum_{r>s} \ln{\textstyle \left(\!\!{n_rn_s \choose e_{rs}}\!\!\right)} + \sum_r \ln{\textstyle \left(\!\!{\left(\!{n_r\choose 2}\!\right)\choose e_{rr}/2}\!\!\right)}\\
            \mathcal{S}^d_m  &= \sum_{rs} \ln{\textstyle \left(\!\!{n_rn_s \choose e_{rs}}\!\!\right)},

479
480
481
        for multigraphs (i.e. ``multigraph == True``). A dense entropy for the
        degree-corrected model is not available, and if requested will return a
        :exc:`NotImplementedError`.
482

483
484
        If ``complete == False`` constants in the above equations which do not
        depend on the partition of the nodes will be omitted.
485

486
487
        Note that in all cases if ``norm==True`` the value returned corresponds
        to the entropy `per edge`, i.e. :math:`(\mathcal{S}_{t/c}\; [\,+\,\mathcal{L}_{t/c}])/ E`.
488
489
        """

490
491
492
        xi_fast = kwargs.get("xi_fast", False)
        dl_deg_alt = kwargs.get("dl_deg_alt", True)

493
494
495
        E = self.E
        N = self.N

496
497
        if dense:
            if self.deg_corr:
498
                raise NotImplementedError('A degree-corrected "dense" entropy is not yet implemented')
499

500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
            S = libcommunity.entropy_dense(self.bg._Graph__graph,
                                           _prop("e", self.bg, self.mrs),
                                           _prop("v", self.bg, self.wr),
                                           multigraph)
        else:
            S = libcommunity.entropy(self.bg._Graph__graph,
                                     _prop("e", self.bg, self.mrs),
                                     _prop("v", self.bg, self.mrp),
                                     _prop("v", self.bg, self.mrm),
                                     _prop("v", self.bg, self.wr),
                                     self.deg_corr)

            if __test__:
                assert not isnan(S) and not isinf(S), "invalid entropy %g (%s) " % (S, str(dict(complete=complete,
                                                                                                random=random, dl=dl,
                                                                                                partition_dl=partition_dl,
                                                                                                dense=dense, multigraph=multigraph,
                                                                                                norm=norm)))
518
519
520
            if complete:
                if self.deg_corr:
                    S -= E
521
522
523
524
                    S += libcommunity.deg_entropy_term(self.g._Graph__graph,
                                                       libcore.any(),
                                                       self.overlap_stats,
                                                       self.N)
525
526
                else:
                    S += E
527

528
529
530
531
532
533
534
535
536
537
                if multigraph:
                    S += libcommunity.entropy_parallel(self.g._Graph__graph,
                                                       _prop("e", self.g, self.eweight))

                if __test__:
                    assert not isnan(S) and not isinf(S), "invalid entropy %g (%s) " % (S, str(dict(complete=complete,
                                                                                                    random=random, dl=dl,
                                                                                                    partition_dl=partition_dl,
                                                                                                    dense=dense, multigraph=multigraph,
                                                                                                    norm=norm)))
538
        if dl:
539
540
541
542
543
544
545
546
547
548
549
550
551
552
            if partition_dl:
                if self.partition_stats.is_enabled():
                    S += self.partition_stats.get_partition_dl()
                else:
                    self.__init_partition_stats(empty=False)
                    S += self.partition_stats.get_partition_dl()
                    self.__init_partition_stats(empty=True)

                if __test__:
                    assert not isnan(S) and not isinf(S), "invalid entropy %g (%s) " % (S, str(dict(complete=complete,
                                                                                                    random=random, dl=dl,
                                                                                                    partition_dl=partition_dl,
                                                                                                    dense=dense, multigraph=multigraph,
                                                                                                    norm=norm)))
553
554
555
            else:
                S += model_entropy(self.B, N, E, directed=self.g.is_directed(), nr=self.wr.a) * E

556
            if self.deg_corr:
557
558
559
560
561
562
                if self.partition_stats.is_enabled():
                    S_seq = self.partition_stats.get_deg_dl(dl_ent, dl_deg_alt, xi_fast)
                else:
                    self.__init_partition_stats(empty=False)
                    S_seq = self.partition_stats.get_deg_dl(dl_ent, dl_deg_alt, xi_fast)
                    self.__init_partition_stats(empty=True)
563

564
                S += S_seq
565

566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
                if __test__:
                    assert not isnan(S_seq) and not isinf(S_seq), "invalid entropy %g (%s) " % (S_seq, str(dict(complete=complete,
                                                                                                                random=random, dl=dl,
                                                                                                                partition_dl=partition_dl,
                                                                                                                dense=dense, multigraph=multigraph,
                                                                                                                norm=norm)))

        if __test__:
            assert not isnan(S) and not isinf(S), "invalid entropy %g (%s) " % (S, str(dict(complete=complete,
                                                                                            random=random, dl=dl,
                                                                                            partition_dl=partition_dl,
                                                                                            dense=dense, multigraph=multigraph,
                                                                                            norm=norm)))

        if norm:
            return S / E
        else:
            return S
584

585
586
587
    def get_matrix(self):
        r"""Returns the block matrix (as a sparse :class:`~scipy.sparse.csr_matrix`),
        which contains the number of edges between each block pair.
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603

        Examples
        --------

        .. testsetup:: get_matrix

           gt.seed_rng(42)
           np.random.seed(42)
           from pylab import *

        .. doctest:: get_matrix

           >>> g = gt.collection.data["polbooks"]
           >>> state = gt.BlockState(g, B=5, deg_corr=True)
           >>> for i in range(1000):
           ...     ds, nmoves = gt.mcmc_sweep(state)
604
           >>> m = state.get_matrix()
605
606
           >>> figure()
           <...>
607
           >>> matshow(m.todense())
608
609
610
611
612
613
614
615
616
617
618
619
620
           <...>
           >>> savefig("bloc_mat.pdf")

        .. testcleanup:: get_matrix

           savefig("bloc_mat.png")

        .. figure:: bloc_mat.*
           :align: center

           A  5x5 block matrix.

       """
621

622
        return adjacency(self.bg, weight=self.mrs)
623
624


625
def model_entropy(B, N, E, directed=False, nr=None):
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
    r"""Computes the amount of information necessary for the parameters of the traditional blockmodel ensemble, for ``B`` blocks, ``N`` vertices, ``E`` edges, and either a directed or undirected graph.

    A traditional blockmodel is defined as a set of :math:`N` vertices which can
    belong to one of :math:`B` blocks, and the matrix :math:`e_{rs}` describes
    the number of edges from block :math:`r` to :math:`s` (or twice that number
    if :math:`r=s` and the graph is undirected).

    For an undirected graph, the number of distinct :math:`e_{rs}` matrices is given by,

    .. math::

       \Omega_m = \left(\!\!{\left(\!{B \choose 2}\!\right) \choose E}\!\!\right)

    and for a directed graph,

    .. math::
       \Omega_m = \left(\!\!{B^2 \choose E}\!\!\right)


    where :math:`\left(\!{n \choose k}\!\right) = {n+k-1\choose k}` is the
    number of :math:`k` combinations with repetitions from a set of size :math:`n`.

    The total information necessary to describe the model is then,

    .. math::

652
653
       \mathcal{L}_t = \ln\Omega_m + \ln\left(\!\!{B \choose N}\!\!\right) + \ln N! - \sum_r \ln n_r!,

654

655
656
    where the remaining term is the information necessary to describe the
    block partition, where :math:`n_r` is the number of nodes in block :math:`r`.
657

658
659
660
661
    If ``nr`` is ``None``, it is assumed :math:`n_r=N/B`.

    The value returned corresponds to the information per edge, i.e.
    :math:`\mathcal{L}_t/E`.
662
663
664
665

    References
    ----------

Tiago Peixoto's avatar
Tiago Peixoto committed
666
667
    .. [peixoto-parsimonious-2013] Tiago P. Peixoto, "Parsimonious module inference in large networks",
       Phys. Rev. Lett. 110, 148701 (2013), :doi:`10.1103/PhysRevLett.110.148701`, :arxiv:`1212.4794`.
Tiago Peixoto's avatar
Tiago Peixoto committed
668
669
670
    .. [peixoto-hierarchical-2014] Tiago P. Peixoto, "Hierarchical block structures and high-resolution
       model selection in large networks ", Phys. Rev. X 4, 011047 (2014), :doi:`10.1103/PhysRevX.4.011047`,
       :arxiv:`1310.4377`.
671
672
673

    """

674
675
676
677
678
679
    if directed:
        x = (B * B);
    else:
        x = (B * (B + 1)) / 2;
    L = lbinom(x + E - 1, E) + partition_entropy(B, N, nr)
    return L / E
680
681
682
683

def Sdl(B, S, N, E, directed=False):
    return S + model_entropy(B, N, E, directed)

684
def lbinom(n, k):
685
    return scipy.special.gammaln(float(n + 1)) - scipy.special.gammaln(float(n - k + 1)) - scipy.special.gammaln(float(k + 1))
686
687
688
689
690

def partition_entropy(B, N, nr=None):
    if nr is None:
        S = N * log(B) + log1p(-(1 - 1./B) ** N)
    else:
691
        S = lbinom(B + N - 1, N) + scipy.special.gammaln(N + 1) - scipy.special.gammaln(nr + 1).sum()
692
    return S
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712

def get_max_B(N, E, directed=False):
    r"""Return the maximum detectable number of blocks, obtained by minimizing:

    .. math::

        \mathcal{L}_t(B, N, E) - E\ln B

    where :math:`\mathcal{L}_t(B, N, E)` is the information necessary to
    describe a traditional blockmodel with `B` blocks, `N` nodes and `E`
    edges (see :func:`model_entropy`).

    Examples
    --------

    >>> gt.get_max_B(N=1e6, E=5e6)
    1572

    References
    ----------
Tiago Peixoto's avatar
Tiago Peixoto committed
713
714
    .. [peixoto-parsimonious-2013] Tiago P. Peixoto, "Parsimonious module inference in large networks",
       Phys. Rev. Lett. 110, 148701 (2013), :doi:`10.1103/PhysRevLett.110.148701`, :arxiv:`1212.4794`.
715
716
717
718
719
720
721
722
723
724
725


    """

    B = fminbound(lambda B: Sdl(B, -log(B), N, E, directed), 1, E,
                  xtol=1e-6, maxfun=1500, disp=0)
    if isnan(B):
        B = 1
    return max(int(ceil(B)), 2)

def get_akc(B, I, N=float("inf"), directed=False):
Tiago Peixoto's avatar
Tiago Peixoto committed
726
    r"""Return the minimum value of the average degree of the network, so that some block structure with :math:`B` blocks can be detected, according to the minimum description length criterion.
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759

    This is obtained by solving

    .. math::

       \Sigma_b = \mathcal{L}_t(B, N, E) - E\mathcal{I}_{t/c} = 0,

    where :math:`\mathcal{L}_{t}` is the necessary information to describe the
    blockmodel parameters (see :func:`model_entropy`), and
    :math:`\mathcal{I}_{t/c}` characterizes the planted block structure, and is
    given by

    .. math::

        \mathcal{I}_t &= \sum_{rs}m_{rs}\ln\left(\frac{m_{rs}}{w_rw_s}\right),\\
        \mathcal{I}_c &= \sum_{rs}m_{rs}\ln\left(\frac{m_{rs}}{m_rm_s}\right),

    where :math:`m_{rs} = e_{rs}/2E` (or :math:`m_{rs} = e_{rs}/E` for directed
    graphs) and :math:`w_r=n_r/N`. We note that :math:`\mathcal{I}_{t/c}\in[0,
    \ln B]`. In the case where :math:`E \gg B^2`, this simplifies to

    .. math::

       \left<k\right>_c &= \frac{2\ln B}{\mathcal{I}_{t/c}},\\
       \left<k^{-/+}\right>_c &= \frac{\ln B}{\mathcal{I}_{t/c}},

    for undirected and directed graphs, respectively. This limit is assumed if
    ``N == inf``.

    Examples
    --------

    >>> gt.get_akc(10, log(10) / 100, N=100)
760
    2.414413200430159
761
762
763

    References
    ----------
Tiago Peixoto's avatar
Tiago Peixoto committed
764
765
    .. [peixoto-parsimonious-2013] Tiago P. Peixoto, "Parsimonious module inference in large networks",
       Phys. Rev. Lett. 110, 148701 (2013), :doi:`10.1103/PhysRevLett.110.148701`, :arxiv:`1212.4794`.
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780

    """
    if N != float("inf"):
        if directed:
            get_dl = lambda ak: model_entropy(B, N, N * ak, directed) - N * ak * I
        else:
            get_dl = lambda ak: model_entropy(B, N, N * ak / 2., directed) - N * ak * I / 2.
        ak = fsolve(lambda ak: get_dl(ak), 10)
        ak = float(ak)
    else:
        ak = 2 * log(B) / S
        if directed:
            ak /= 2
    return ak

781
782
783
784
def mcmc_sweep(state, beta=1., c=1., dl=False, dense=False, multigraph=False,
               node_coherent=False, nmerges=0, nmerge_sweeps=1, merge_map=None,
               coherent_merge=False, sequential=True, parallel=False,
               vertices=None, verbose=False, **kwargs):
785
    r"""Performs a Markov chain Monte Carlo 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.
786
787
788

    Parameters
    ----------
789
    state : :class:`~graph_tool.community.BlockState` or :class:`~graph_tool.community.OverlapBlockState`
790
        The block state.
791
    beta : ``float`` (optional, default: `1.0`)
792
        The inverse temperature parameter :math:`\beta`.
793
794
795
796
797
    c : ``float`` (optional, default: ``1.0``)
        This parameter specifies how often fully random moves are attempted,
        instead of more likely moves based on the inferred block partition.
        For ``c == 0``, no fully random moves are attempted, and for ``c == inf``
        they are always attempted.
798
799
800
    dl : ``bool`` (optional, default: ``False``)
        If ``True``, the change in the whole description length will be
        considered after each vertex move, not only the entropy.
801
802
803
804
805
    dense : ``bool`` (optional, default: ``False``)
        If ``True``, the "dense" variant of the entropy will be computed.
    multigraph : ``bool`` (optional, default: ``False``)
        If ``True``, the multigraph entropy will be used. Only has an effect
        if ``dense == True``.
806
807
808
809
    node_coherent : ``bool`` (optional, default: ``False``)
        If ``True``, and if the ``state`` is an instance of
        :class:`~graph_tool.community.OverlapBlockState`, then all half-edges
        incident on the same node are moved simultaneously.
810
811
812
813
814
    sequential : ``bool`` (optional, default: ``True``)
        If ``True``, the move attempts on the vertices are done in sequential
        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.
815
816
817
818
819
820
821
822
823
824
    parallel : ``bool`` (optional, default: ``False``)
        If ``True``, the updates are performed in parallel (multiple
        threads).

        .. warning::

            If this is used, the Markov Chain is not guaranteed to be sampled with
            the correct probabilities. This is better used in conjunction with
            ``beta=float('inf')``, where this is not an issue.

Tiago Peixoto's avatar
Tiago Peixoto committed
825
    vertices : ``list of ints`` (optional, default: ``None``)
826
827
        A list of vertices which will be attempted to be moved. If ``None``, all
        vertices will be attempted.
828
829
830
831
832
833
    verbose : ``bool`` (optional, default: ``False``)
        If ``True``, verbose information is displayed.

    Returns
    -------

834
    dS : ``float``
835
836
837
838
839
840
841
842
       The entropy difference (per edge) after a full sweep.
    nmoves : ``int``
       The number of accepted block membership moves.


    Notes
    -----

843
    This algorithm performs a Markov chain Monte Carlo sweep on the network,
844
845
    where the block memberships are randomly moved, and either accepted or
    rejected, so that after sufficiently many sweeps the partitions are sampled
846
    with probability proportional to :math:`e^{-\beta\mathcal{S}_{t/c}}`, where
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
    :math:`\mathcal{S}_{t/c}` is the blockmodel entropy, given by

    .. math::

      \mathcal{S}_t &\cong - \frac{1}{2} \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{n_rn_s}\right), \\
      \mathcal{S}^d_t &\cong - \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{n_rn_s}\right),

    for undirected and directed traditional blockmodels (``deg_corr == False``),
    respectively, where :math:`e_{rs}` is the number of edges from block
    :math:`r` to :math:`s` (or the number of half-edges for the undirected case
    when :math:`r=s`), and :math:`n_r` is the number of vertices in block
    :math:`r`, and constant terms which are independent of the block partition
    were dropped (see :meth:`BlockState.entropy` for the complete entropy). For
    the degree-corrected variant with "hard" degree constraints the equivalent
    expressions are

    .. math::

       \mathcal{S}_c &\cong  - \frac{1}{2} \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{e_re_s}\right), \\
       \mathcal{S}^d_c &\cong - \sum_{rs}e_{rs}\ln\left(\frac{e_{rs}}{e^+_re^-_s}\right),

    where :math:`e_r = \sum_se_{rs}` is the number of half-edges incident on
    block :math:`r`, and :math:`e^+_r = \sum_se_{rs}` and :math:`e^-_r =
    \sum_se_{sr}` are the number of out- and in-edges adjacent to block
    :math:`r`, respectively.

    The Monte Carlo algorithm employed attempts to improve the mixing time of
874
875
    the Markov chain by proposing membership moves :math:`r\to s` with
    probability :math:`p(r\to s|t) \propto e_{ts} + c`, where :math:`t` is the
876
    block label of a random neighbour of the vertex being moved. See
877
    [peixoto-efficient-2014]_ for more details.
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918

    This algorithm has a complexity of :math:`O(E)`, where :math:`E` is the
    number of edges in the network.

    Examples
    --------
    .. testsetup:: mcmc

       gt.seed_rng(42)
       np.random.seed(42)

    .. doctest:: mcmc

       >>> g = gt.collection.data["polbooks"]
       >>> state = gt.BlockState(g, B=3, deg_corr=True)
       >>> pv = None
       >>> for i in range(1000):        # remove part of the transient
       ...     ds, nmoves = gt.mcmc_sweep(state)
       >>> for i in range(1000):
       ...     ds, nmoves = gt.mcmc_sweep(state)
       ...     pv = gt.collect_vertex_marginals(state, pv)
       >>> gt.graph_draw(g, pos=g.vp["pos"], vertex_shape="pie", vertex_pie_fractions=pv, output="polbooks_blocks_soft.pdf")
       <...>

    .. testcleanup:: mcmc

       gt.graph_draw(g, pos=g.vp["pos"], vertex_shape="pie", vertex_pie_fractions=pv, output="polbooks_blocks_soft.png")

    .. figure:: polbooks_blocks_soft.*
       :align: center

       "Soft" block partition of a political books network with :math:`B=3`.

     References
    ----------

    .. [holland-stochastic-1983] Paul W. Holland, Kathryn Blackmond Laskey,
       Samuel Leinhardt, "Stochastic blockmodels: First steps",
       Carnegie-Mellon University, Pittsburgh, PA 15213, U.S.A., :doi:`10.1016/0378-8733(83)90021-7`
    .. [faust-blockmodels-1992] Katherine Faust, and Stanley
       Wasserman. "Blockmodels: Interpretation and Evaluation." Social Networks
919
       14, no. 1-2 (1992): 5-61. :doi:`10.1016/0378-8733(92)90013-W`
920
921
922
923
924
925
    .. [karrer-stochastic-2011] Brian Karrer, and M. E. J. Newman. "Stochastic
       Blockmodels and Community Structure in Networks." Physical Review E 83,
       no. 1 (2011): 016107. :doi:`10.1103/PhysRevE.83.016107`.
    .. [peixoto-entropy-2012] Tiago P. Peixoto "Entropy of Stochastic Blockmodel
       Ensembles." Physical Review E 85, no. 5 (2012): 056122. :doi:`10.1103/PhysRevE.85.056122`,
       :arxiv:`1112.6028`.
Tiago Peixoto's avatar
Tiago Peixoto committed
926
927
    .. [peixoto-parsimonious-2013] Tiago P. Peixoto, "Parsimonious module inference in large networks",
       Phys. Rev. Lett. 110, 148701 (2013), :doi:`10.1103/PhysRevLett.110.148701`, :arxiv:`1212.4794`.
928
929
930
    .. [peixoto-efficient-2014] Tiago P. Peixoto, "Efficient Monte Carlo and greedy
       heuristic for the inference of stochastic block models", Phys. Rev. E 89, 012804 (2014),
       :doi:`10.1103/PhysRevE.89.012804`, :arxiv:`1310.4378`.
931
932
933
    .. [peixoto-model-2014] Tiago P. Peixoto, "Model selection and hypothesis
       testing for large-scale network models with overlapping groups",
       :arxiv:`1409.3059`.
934
935
    """

936
    if state.B == 1:
937
938
        return 0., 0

939
    if vertices is not None:
940
941
942
        vlist = libcommunity.get_vector(len(vertices))
        vlist.a = vertices
        vertices = vlist
943
        #state.sweep_vertices = vertices
944

945
946
947
948
949
    if state.sweep_vertices is None:
        vertices = libcommunity.get_vector(state.g.num_vertices())
        vertices.a = state.g.vertex_index.copy("int").fa
        state.sweep_vertices = vertices

950
951
    random_move = c == float("inf")

952
    if random_move or nmerges > 0:
953
954
955
956
        state._BlockState__build_egroups(empty=True)
    elif state.egroups is None:
        state._BlockState__build_egroups(empty=False)

957
    bclabel = state.get_bclabel()
958

959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
    if nmerges == 0:
        nmerge_sweeps = 1
        if state.nsampler is None:
            state._BlockState__build_nsampler(empty=state.overlap)
        nsampler = state.nsampler
        ncavity_sampler = state.nsampler
        merge_map = state.g.new_vertex_property("int")
    else:
        if kwargs.get("unweighted_merge", False):
            emask = state.mrs
        else:
            emask = state.mrs.copy()
            emask.a = emask.a > 0

        nsampler = libcommunity.init_neighbour_sampler(state.bg._Graph__graph,
                                                       _prop("e", state.bg, emask),
                                                       True, False)
        ncavity_sampler = libcommunity.init_neighbour_sampler(state.bg._Graph__graph,
                                                              _prop("e", state.bg, emask),
                                                              False, False)
        beta = float("inf")

        if merge_map is None:
            merge_map = state.g.vertex_index.copy("int")

    if dl and not state.partition_stats.is_enabled():
        if state.overlap:
            state._OverlapBlockState__init_partition_stats(empty=False)
        else:
            state._BlockState__init_partition_stats(empty=False)

    if not dl and state.partition_stats.is_enabled():
        if state.overlap:
            state._OverlapBlockState__init_partition_stats(empty=True)
        else:
            state._BlockState__init_partition_stats(empty=True)

    if __test__:
        assert state._BlockState__check_clabel(), "clabel already invalid!"
        S = state.entropy(dense=dense, multigraph=multigraph, complete=False, dl=dl, dl_deg_alt=False, xi_fast=True)
        assert not (isinf(S) or isnan(S)), "invalid entropy before sweep: %g" % S
1000

For faster browsing, not all history is shown. View entire blame