base_states.py 29.3 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-2022 Tiago de Paula Peixoto <tiago@skewed.de>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU Lesser 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

21
from .. import Vector_size_t, group_vector_property
22
23
24
25
26
27
28
29
30

from . util import *

import functools
from abc import ABC, abstractmethod
import gc

import numpy

31
32
import graph_tool.draw

33
34
35
36
37
38
39
40
41
42
43
44
__test__ = False

def set_test(test):
    global __test__
    __test__ = test

def _bm_test():
    global __test__
    return __test__

def copy_state_wrap(func):
    @functools.wraps(func)
45
    def wrapper(self, *args, test=True, **kwargs):
46
47
48

        S = func(self, *args, **kwargs)

49
        if _bm_test() and test:
50
51
52
53
            assert not isnan(S) and not isinf(S), \
                "invalid entropy %g (%s) " % (S, str(args))

            state_copy = self.copy()
54
            Salt = state_copy.entropy(*args, test=False, **kwargs)
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

            assert math.isclose(S, Salt, abs_tol=1e-8), \
                "entropy discrepancy after copying (%g %g %g)" % (S, Salt,
                                                                  S - Salt)
        return S

    return wrapper


def mcmc_sweep_wrap(func):
    @functools.wraps(func)
    def wrapper(self, *args, **kwargs):
        test = kwargs.pop("test", True)
        entropy_args = kwargs.get("entropy_args", {})

70
71
72
        if not kwargs.get("dispatch", True):
            return func(self, *args, **kwargs)

73
74
75
76
77
        if _bm_test() and test:
            if hasattr(self, "_check_clabel"):
                assert self._check_clabel(), "invalid clabel before sweep"
            Si = self.entropy(**entropy_args)

78
        ret = func(self, *args, **kwargs)
79
80
81
82

        if _bm_test() and test:
            if hasattr(self, "_check_clabel"):
                assert self._check_clabel(), "invalid clabel after sweep"
83
            dS = ret[0]
84
85
86
87
            Sf = self.entropy(**entropy_args)
            assert math.isclose(dS, (Sf - Si), abs_tol=1e-8), \
                "inconsistent entropy delta %g (%g): %s" % (dS, Sf - Si,
                                                            str(entropy_args))
88
        return ret
89
90
91
    return wrapper

class MCMCState(ABC):
Tiago Peixoto's avatar
Tiago Peixoto committed
92
    r"""Base state that implements single-flip MCMC sweeps"""
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

    @abstractmethod
    def _mcmc_sweep_dispatch(self, mcmc_state):
        pass

    @abstractmethod
    def _get_entropy_args(self, kwargs):
        pass

    @mcmc_sweep_wrap
    def mcmc_sweep(self, beta=1., c=.5, d=.01, niter=1, entropy_args={},
                   allow_vacate=True, sequential=True, deterministic=False,
                   vertices=None, verbose=False, **kwargs):
        r"""Perform ``niter`` sweeps of a Metropolis-Hastings acceptance-rejection
        MCMC to sample network partitions.

        Parameters
        ----------
        beta : ``float`` (optional, default: ``1.``)
            Inverse temperature.
        c : ``float`` (optional, default: ``.5``)
            Sampling parameter ``c`` for move proposals: For :math:`c\to 0` the
            blocks are sampled according to the local neighborhood of a given
            node and their block connections; for :math:`c\to\infty` the blocks
            are sampled randomly. Note that only for :math:`c > 0` the MCMC is
            guaranteed to be ergodic.
        d : ``float`` (optional, default: ``.01``)
            Probability of selecting a new (i.e. empty) group for a given move.
        niter : ``int`` (optional, default: ``1``)
            Number of sweeps to perform. During each sweep, a move attempt is
            made for each node.
        entropy_args : ``dict`` (optional, default: ``{}``)
            Entropy arguments, with the same meaning and defaults as in
126
            :meth:`graph_tool.inference.BlockState.entropy`.
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
        allow_vacate : ``bool`` (optional, default: ``True``)
            Allow groups to be vacated.
        sequential : ``bool`` (optional, default: ``True``)
            If ``sequential == True`` each vertex move attempt is made
            sequentially, where vertices are visited in random order. Otherwise
            the moves are attempted by sampling vertices randomly, so that the
            same vertex can be moved more than once, before other vertices had
            the chance to move.
        deterministic : ``bool`` (optional, default: ``False``)
            If ``sequential == True`` and ``deterministic == True`` the
            vertices will be visited in deterministic order.
        vertices : ``list`` of ints (optional, default: ``None``)
            If provided, this should be a list of vertices which will be
            moved. Otherwise, all vertices will.
        verbose : ``bool`` (optional, default: ``False``)
            If ``verbose == True``, detailed information will be displayed.

        Returns
        -------
        dS : ``float``
            Entropy difference after the sweeps.
        nattempts : ``int``
            Number of vertex moves attempted.
        nmoves : ``int``
            Number of vertices moved.

        Notes
        -----
        This algorithm has an :math:`O(E)` complexity, where :math:`E` is the
Tiago Peixoto's avatar
Tiago Peixoto committed
156
        number of edges (independent of the number of groups).
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189

        References
        ----------
        .. [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`
        """

        mcmc_state = DictState(locals())
        mcmc_state.oentropy_args = self._get_entropy_args(entropy_args)
        mcmc_state.vlist = Vector_size_t()
        if vertices is None:
            vertices = self.g.vertex_index.copy().fa
            if getattr(self, "is_weighted", False):
                # ignore vertices with zero weight
                vw = self.vweight.fa
                vertices = vertices[vw > 0]
        mcmc_state.vlist.resize(len(vertices))
        mcmc_state.vlist.a = vertices
        mcmc_state.state = self._state

        dispatch = kwargs.pop("dispatch", True)

        if len(kwargs) > 0:
            raise ValueError("unrecognized keyword arguments: " +
                             str(list(kwargs.keys())))
        if dispatch:
            return self._mcmc_sweep_dispatch(mcmc_state)
        else:
            return mcmc_state

class MultiflipMCMCState(ABC):
Tiago Peixoto's avatar
Tiago Peixoto committed
190
    r"""Base state that implements multiflip (merge-split) MCMC sweeps"""
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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

    @abstractmethod
    def _multiflip_mcmc_sweep_dispatch(self, mcmc_state):
        pass

    @abstractmethod
    def _get_entropy_args(self, kwargs):
        pass

    @mcmc_sweep_wrap
    def multiflip_mcmc_sweep(self, beta=1., c=.5, psingle=None, psplit=1,
                             pmerge=1, pmergesplit=1, d=0.01, gibbs_sweeps=10,
                             niter=1, entropy_args={}, accept_stats=None,
                             verbose=False, **kwargs):
        r"""Perform ``niter`` sweeps of a Metropolis-Hastings acceptance-rejection MCMC
        with multiple simultaneous moves (i.e. merges and splits) to sample
        network partitions.

        Parameters
        ----------
        beta : ``float`` (optional, default: ``1.``)
            Inverse temperature.
        c : ``float`` (optional, default: ``.5``)
            Sampling parameter ``c`` for move proposals: For :math:`c\to 0` the
            blocks are sampled according to the local neighborhood of a given
            node and their block connections; for :math:`c\to\infty` the blocks
            are sampled randomly. Note that only for :math:`c > 0` the MCMC is
            guaranteed to be ergodic.
        psingle : ``float`` (optional, default: ``None``)
            Relative probability of proposing a single node move. If ``None``,
            it will be selected as the number of nodes in the graph.
        psplit : ``float`` (optional, default: ``1``)
            Relative probability of proposing a group split.
        pmergesplit : ``float`` (optional, default: ``1``)
            Relative probability of proposing a marge-split move.
        d : ``float`` (optional, default: ``1``)
            Probability of selecting a new (i.e. empty) group for a given
            single-node move.
        gibbs_sweeps : ``int`` (optional, default: ``10``)
            Number of sweeps of Gibbs sampling to be performed (i.e. each node
            is attempted once per sweep) to refine a split proposal.
        niter : ``int`` (optional, default: ``1``)
            Number of sweeps to perform. During each sweep, a move attempt is
            made for each node, on average.
        entropy_args : ``dict`` (optional, default: ``{}``)
            Entropy arguments, with the same meaning and defaults as in
237
            :meth:`graph_tool.inference.BlockState.entropy`.
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
        verbose : ``bool`` (optional, default: ``False``)
            If ``verbose == True``, detailed information will be displayed.

        Returns
        -------
        dS : ``float``
            Entropy difference after the sweeps.
        nattempts : ``int``
            Number of vertex moves attempted.
        nmoves : ``int``
            Number of vertices moved.

        Notes
        -----
        This algorithm has an :math:`O(E)` complexity, where :math:`E` is the
        number of edges (independent of the number of groups).

        References
        ----------
        .. [peixoto-merge-split-2020] Tiago P. Peixoto, "Merge-split Markov
           chain Monte Carlo for community detection", Phys. Rev. E 102, 012305
           (2020), :doi:`10.1103/PhysRevE.102.012305`, :arxiv:`2003.07070`

        """

        if psingle is None:
            psingle = self.g.num_vertices()
        gibbs_sweeps = max(gibbs_sweeps, 1)
        nproposal = Vector_size_t(4)
        nacceptance = Vector_size_t(4)
        force_move = kwargs.pop("force_move", False)
        mcmc_state = DictState(locals())
        mcmc_state.oentropy_args = self._get_entropy_args(entropy_args)
        mcmc_state.state = self._state

        dispatch = kwargs.pop("dispatch", True)

        if len(kwargs) > 0:
            raise ValueError("unrecognized keyword arguments: " +
                             str(list(kwargs.keys())))
        if dispatch:
            dS, nattempts, nmoves = self._multiflip_mcmc_sweep_dispatch(mcmc_state)
        else:
            return mcmc_state

        if accept_stats is not None:
            for key in ["proposal", "acceptance"]:
                if key not in accept_stats:
                    accept_stats[key] = numpy.zeros(len(nproposal),
                                                    dtype="uint64")
            accept_stats["proposal"] += nproposal.a
            accept_stats["acceptance"] += nacceptance.a

        return dS, nattempts, nmoves

class MultilevelMCMCState(ABC):
Tiago Peixoto's avatar
Tiago Peixoto committed
294
    r"""Base state that implements multilevel agglomerative MCMC sweeps"""
295
296
297
298
299
300
301
302
303

    @abstractmethod
    def _multilevel_mcmc_sweep_dispatch(self, mcmc_state):
        pass

    @abstractmethod
    def _get_entropy_args(self, kwargs):
        pass

304
    def _get_bclabel(self):
305
306
        return None

307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
    @mcmc_sweep_wrap
    def multilevel_mcmc_sweep(self, niter=1, beta=1., c=.5, psingle=None,
                              pmultilevel=1, d=0.01, r=0.9, random_bisect=True,
                              merge_sweeps=10, mh_sweeps=10, init_r=0.99,
                              init_beta=1., gibbs=False, B_min=1,
                              B_max=numpy.iinfo(numpy.uint64).max, b_min=None,
                              b_max=None, M=None, cache_states=True,
                              entropy_args={}, verbose=False, **kwargs):
        r"""Perform ``niter`` sweeps of a multilevel agglomerative acceptance-rejection
        MCMC to sample network partitions, that uses a bisection search on the
        number of groups, together with group merges and singe-node moves.

        Parameters
        ----------
        niter : ``int`` (optional, default: ``1``)
            Number of sweeps to perform. During each sweep, a move attempt is
            made for each node, on average.
        beta : ``float`` (optional, default: ``1.``)
            Inverse temperature.
        c : ``float`` (optional, default: ``.5``)
            Sampling parameter ``c`` for move proposals: For :math:`c\to 0` the
            blocks are sampled according to the local neighborhood of a given
            node and their block connections; for :math:`c\to\infty` the blocks
            are sampled randomly. Note that only for :math:`c > 0` the MCMC is
            guaranteed to be ergodic.
        psingle : ``float`` (optional, default: ``None``)
            Relative probability of proposing a single node move. If ``None``,
            it will be selected as the number of nodes in the graph.
        pmultilevel : ``float`` (optional, default: ``1``)
            Relative probability of proposing a multilevel move.
337
        d : ``float`` (optional, default: ``.01``)
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
            Probability of selecting a new (i.e. empty) group for a given
            single-node move.
        r : ``float`` (optional, default: ``0.9``)
            Group shrink ratio. The number of groups is reduced by this fraction
            at each merge sweep.
        random_bisect : ``bool`` (optional, default: ``True``)
            If ``True``, bisections are done at randomly chosen
            intervals. Otherwise a Fibonacci sequence is used.
        merge_sweeps : ``int`` (optional, default: ``10``)
            Number of sweeps spent to find good merge proposals.
        mh_sweeps : ``int`` (optional, default: ``10``)
            Number of single-node Metropolis-Hastings sweeps between merge splits.
        init_r : ``double`` (optional, default: ``0.99``)
            Stopping criterion for the intialization phase, after each node is
            put in their own group, to set the initial upper bound of the
            bisection search. A number of single-node Metropolis-Hastings sweeps
            is done until the number of groups is shrunk by a factor that is
            larger than this parameter.
        init_beta : ``float`` (optional, default: ``1.``)
            Inverse temperature to be used for the very first sweep of the
            initialization phase.
        gibbs : ``bool`` (optional, default: ``False``)
            If ``True``, the single node moves use (slower) Gibbs sampling,
            rather than Metropolis-Hastings.
        B_min : ``int`` (optional, default: ``1``)
            Minimum number of groups to be considered in the search.
        b_min : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
            If provided, this will be used for the partition corresponding to ``B_min``.
        B_max : ``int`` (optional, default: ``1``)
            Maximum number of groups to be considered in the search.
        b_max : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
            If provided, this will be used for the partition corresponding to ``B_max``.
        M : ``int`` (optional, default: ``None``)
            Maximum number of groups to select for the multilevel move. If
            ``None`` is provided, then all groups are always elected.
        cache_states : ``bool`` (optional, default: ``True``)
            If ``True``, intermediary states will be cached during the bisection
            search.
        entropy_args : ``dict`` (optional, default: ``{}``)
            Entropy arguments, with the same meaning and defaults as in
378
            :meth:`graph_tool.inference.BlockState.entropy`.
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
        verbose : ``bool`` (optional, default: ``False``)
            If ``verbose == True``, detailed information will be displayed.

        Returns
        -------
        dS : ``float``
            Entropy difference after the sweeps.
        nattempts : ``int``
            Number of vertex moves attempted.
        nmoves : ``int``
            Number of vertices moved.

        Notes
        -----
        This algorithm has an :math:`O(E\ln^2 N)` complexity, where :math:`E` is
        the number of edges and :math:`N` is the number of nodes (independently of
        the number of groups).

        References
        ----------
        .. [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`

        """

        if psingle is None:
            psingle = self.g.num_vertices()
        merge_sweeps = max(merge_sweeps, 1)
        if M is None:
            M = self.g.num_vertices()
            global_moves = True
        else:
            global_moves = False
414
415
416
        bclabel = self._get_bclabel()
        if bclabel is not None:
            B_min = max(len(numpy.unique(bclabel.fa)), B_min)
417
        if b_min is None:
418
            b_min = self.g.vertex_index.copy("int")
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
        if b_max is None:
            b_max = self.g.new_vp("int")

        mcmc_state = DictState(locals())
        mcmc_state.oentropy_args = self._get_entropy_args(entropy_args)
        mcmc_state.state = self._state

        dispatch = kwargs.pop("dispatch", True)

        if len(kwargs) > 0:
            raise ValueError("unrecognized keyword arguments: " +
                             str(list(kwargs.keys())))
        if dispatch:
            return self._multilevel_mcmc_sweep_dispatch(mcmc_state)
        else:
            return mcmc_state

class GibbsMCMCState(ABC):
Tiago Peixoto's avatar
Tiago Peixoto committed
437
    r"""Base state that implements single flip MCMC sweeps"""
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462

    @abstractmethod
    def _gibbs_sweep_dispatch(self, gibbs_state):
        pass

    @abstractmethod
    def _get_entropy_args(self, kwargs):
        pass

    @mcmc_sweep_wrap
    def gibbs_sweep(self, beta=1., niter=1, entropy_args={},
                    allow_new_group=True, sequential=True, deterministic=False,
                    vertices=None, verbose=False, **kwargs):
        r"""Perform ``niter`` sweeps of a rejection-free Gibbs MCMC to sample network
        partitions.

        Parameters
        ----------
        beta : ``float`` (optional, default: ``1.``)
            Inverse temperature.
        niter : ``int`` (optional, default: ``1``)
            Number of sweeps to perform. During each sweep, a move attempt is
            made for each node.
        entropy_args : ``dict`` (optional, default: ``{}``)
            Entropy arguments, with the same meaning and defaults as in
463
            :meth:`graph_tool.inference.BlockState.entropy`.
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
        allow_new_group : ``bool`` (optional, default: ``True``)
            Allow the number of groups to increase and decrease.
        sequential : ``bool`` (optional, default: ``True``)
            If ``sequential == True`` each vertex move attempt is made
            sequentially, where vertices are visited in random order. Otherwise
            the moves are attempted by sampling vertices randomly, so that the
            same vertex can be moved more than once, before other vertices had
            the chance to move.
        deterministic : ``bool`` (optional, default: ``False``)
            If ``sequential == True`` and ``deterministic == True`` the
            vertices will be visited in deterministic order.
        vertices : ``list`` of ints (optional, default: ``None``)
            If provided, this should be a list of vertices which will be
            moved. Otherwise, all vertices will.
        verbose : ``bool`` (optional, default: ``False``)
            If ``verbose == True``, detailed information will be displayed.

        Returns
        -------
        dS : ``float``
            Entropy difference after the sweeps.
        nattempts : ``int``
            Number of vertex moves attempted.
        nmoves : ``int``
            Number of vertices moved.

        Notes
        -----
        This algorithm has an :math:`O(E\times B)` complexity, where :math:`B`
Tiago Peixoto's avatar
Tiago Peixoto committed
493
        is the number of groups, and :math:`E` is the number of edges.
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514

        """

        gibbs_state = DictState(locals())
        gibbs_state.oentropy_args = self._get_entropy_args(entropy_args)
        gibbs_state.vlist = Vector_size_t()
        if vertices is None:
            vertices = self.g.get_vertices()
        gibbs_state.vlist.resize(len(vertices))
        gibbs_state.vlist.a = vertices
        gibbs_state.state = self._state

        if len(kwargs) > 0:
            raise ValueError("unrecognized keyword arguments: " +
                             str(list(kwargs.keys())))

        dS, nattempts, nmoves = self._gibbs_sweep_dispatch(gibbs_state)

        return dS, nattempts, nmoves

class MulticanonicalMCMCState(ABC):
Tiago Peixoto's avatar
Tiago Peixoto committed
515
    r"""Base state that implements multicanonical MCMC sweeps"""
516
517
518
519
520
521
522
523
524
525
526
527

    @abstractmethod
    def _multicanonical_sweep_dispatch(self, multicanonical_state):
        pass

    @mcmc_sweep_wrap
    def multicanonical_sweep(self, m_state, multiflip=False, **kwargs):
        r"""Perform ``niter`` sweeps of a non-Markovian multicanonical sampling using the
        Wang-Landau algorithm.

        Parameters
        ----------
528
529
        m_state : :class:`~graph_tool.inference.MulticanonicalState`
            :class:`~graph_tool.inference.MulticanonicalState` instance
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
            containing the current state of the Wang-Landau run.
        multiflip : ``bool`` (optional, default: ``False``)
            If ``True``, ``multiflip_mcmc_sweep()`` will be used, otherwise
            ``mcmc_sweep()``.
        **kwargs : Keyword parameter list
            The remaining parameters will be passed to
            ``multiflip_mcmc_sweep()`` or ``mcmc_sweep()``.

        Returns
        -------
        dS : ``float``
            Entropy difference after the sweeps.
        nattempts : ``int``
            Number of vertex moves attempted.
        nmoves : ``int``
            Number of vertices moved.

        Notes
        -----
        This algorithm has an :math:`O(E)` complexity, where :math:`E` is the
Tiago Peixoto's avatar
Tiago Peixoto committed
550
        number of edges (independent of the number of groups).
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567

        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`
        """

        if not multiflip:
            kwargs["sequential"] = False
        kwargs["beta"] = 1

        args = dmask(locals(), ["self", "kwargs"])
        multi_state = DictState(args)

        entropy_args = kwargs.get("entropy_args", {})
568
        entropy_offset = kwargs.pop("entropy_offset", 0)
569
570
571
572
573
574
575
576
577

        if multiflip:
            mcmc_state = self.multiflip_mcmc_sweep(dispatch=False, **kwargs)
        else:
            mcmc_state = self.mcmc_sweep(dispatch=False, **kwargs)

        multi_state.update(mcmc_state)
        multi_state.multiflip = multiflip

578
        multi_state.S = self.entropy(**entropy_args) + entropy_offset
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
        multi_state.state = self._state

        multi_state.f = m_state._f
        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

        if (multi_state.S < multi_state.S_min or
            multi_state.S > multi_state.S_max):
            raise ValueError("initial entropy %g out of bounds (%g, %g)" %
                             (multi_state.S, multi_state.S_min,
                              multi_state.S_max))

        S, nattempts, nmoves = self._multicanonical_sweep_dispatch(multi_state)

        return S, nattempts, nmoves

Tiago Peixoto's avatar
Tiago Peixoto committed
597
598
class ExhaustiveSweepState(ABC):
    r"""Base state that implements exhaustive enumerative sweeps"""
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615

    @abstractmethod
    def _exhaustive_sweep_dispatch(self, exhaustive_state):
        pass

    @abstractmethod
    def _get_entropy_args(self, kwargs):
        pass

    def exhaustive_sweep(self, entropy_args={}, callback=None, density=None,
                         vertices=None, initial_partition=None, max_iter=None):
        r"""Perform an exhaustive loop over all possible network partitions.

        Parameters
        ----------
        entropy_args : ``dict`` (optional, default: ``{}``)
            Entropy arguments, with the same meaning and defaults as in
616
            :meth:`graph_tool.inference.BlockState.entropy`.
617
618
619
620
621
622
623
624
625
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
652
653
654
655
        callback : callable object (optional, default: ``None``)
            Function to be called for each partition, with three arguments ``(S,
            S_min, b_min)`` corresponding to the the current entropy value, the
            minimum entropy value so far, and the corresponding partition,
            respectively. If not provided, and ``hist is None`` an iterator over
            the same values will be returned instead.
        density : ``tuple`` (optional, default: ``None``)
            If provided, it should contain a tuple with values ``(S_min, S_max,
            n_bins)``, which will be used to obtain the density of states via a
            histogram of size ``n_bins``. This parameter is ignored unless
            ``callback is None``.
        vertices : iterable of ints (optional, default: ``None``)
            If provided, this should be a list of vertices which will be
            moved. Otherwise, all vertices will.
        initial_partition : iterable  of ints (optional, default: ``None``)
            If provided, this will provide the initial partition for the
            iteration.
        max_iter : ``int`` (optional, default: ``None``)
            If provided, this will limit the total number of iterations.

        Returns
        -------
        states : iterator over (S, S_min, b_min)
            If ``callback`` is ``None`` and ``hist`` is ``None``, the function
            will return an iterator over ``(S, S_min, b_min)`` corresponding to
            the the current entropy value, the minimum entropy value so far, and
            the corresponding partition, respectively.
        Ss, counts : pair of :class:`numpy.ndarray`
            If ``callback is None`` and ``hist is not None``, the function will
            return the values of each bin (``Ss``) and the state count of each
            bin (``counts``).
        b_min : :class:`~graph_tool.VertexPropertyMap`
            If ``callback is not None`` or ``hist is not None``, the function
            will also return partition with smallest entropy.

        Notes
        -----

        This algorithm has an :math:`O(B^N)` complexity, where :math:`B` is the
Tiago Peixoto's avatar
Tiago Peixoto committed
656
        number of groups, and :math:`N` is the number of vertices.
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694

        """

        exhaustive_state = DictState(dict(max_iter=max_iter if max_iter is not None else 0))
        exhaustive_state.oentropy_args = self._get_entropy_args(entropy_args)
        exhaustive_state.vlist = Vector_size_t()
        if vertices is None:
            vertices = self.g.vertex_index.copy().fa
            if getattr(self, "is_weighted", False):
                # ignore vertices with zero weight
                vw = self.vweight.fa
                vertices = vertices[vw > 0]
        if initial_partition is None:
            initial_partition = zeros(len(vertices), dtype="uint64")
        self.move_vertex(vertices, initial_partition)
        exhaustive_state.vlist.resize(len(vertices))
        exhaustive_state.vlist.a = vertices
        exhaustive_state.S = self.entropy(**entropy_args)
        exhaustive_state.state = self._state
        exhaustive_state.b_min = b_min = self.g.new_vp("int32_t")

        if density is not None:
            density = (density[0], density[1],
                       numpy.zeros(density[2], dtype="uint64"))
        if callback is not None:
            _callback = lambda S, S_min: callback(S, S_min, b_min)
        else:
            _callback = None
        ret = self._exhaustive_sweep_dispatch(exhaustive_state, _callback,
                                              density)
        if _callback is None:
            if density is None:
                return ((S, S_min, b_min) for S, S_min in ret)
            else:
                Ss = numpy.linspace(density[0], density[1], len(density[2]))
                return (Ss, density[2]), b_min
        else:
            return b_min
695
696
697
698
699
700
701
702
703
704


class DrawBlockState(ABC):
    r"""Base state that implements group-based drawing."""

    def draw(self, **kwargs):
        r"""Convenience wrapper to :func:`~graph_tool.draw.graph_draw` that
        draws the state of the graph as colors on the vertices and edges."""
        gradient = self.g.new_ep("double")
        gradient = group_vector_property([gradient])
705
        b = self.g.own_property(self.b)
706
        return graph_tool.draw.graph_draw(self.g,
707
708
709
                                          vertex_fill_color=kwargs.get("vertex_fill_color", b),
                                          vertex_color=kwargs.get("vertex_color", b),
                                          edge_gradient=kwargs.get("edge_gradient", gradient),
710
711
712
                                          **dmask(kwargs, ["vertex_fill_color",
                                                           "vertex_color",
                                                           "edge_gradient"]))