minimize.py 8.29 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
# 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.
12
#
13
14
15
16
# 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.
17
#
18
19
# 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/>.
20
21
22
23
24
25

import numpy
from . util import *
from . blockmodel import *
from . nested_blockmodel import *

26
def minimize_blockmodel_dl(g, state=BlockState, state_args={},
27
                           multilevel_mcmc_args={}):
28
    r"""Fit the stochastic block model, by minimizing its description length using an
29
    agglomerative heuristic.
30
31
32
33
34

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        The graph.
35
36
    state : SBM-like state class (optional, default: :class:`~graph_tool.inference.BlockState`)
        Type of model that will be used. Must be derived from :class:`~graph_tool.inference.MultilevelMCMCState`.
37
38
    state_args : ``dict`` (optional, default: ``{}``)
        Arguments to be passed to appropriate state constructor (e.g.
39
        :class:`~graph_tool.inference.BlockState`)
40
    multilevel_mcmc_args : ``dict`` (optional, default: ``{}``)
41
        Arguments to be passed to :meth:`~graph_tool.inference.MultilevelMCMCState.multilevel_mcmc_sweep`.
42
43
44

    Returns
    -------
45
46
    min_state : type given by parameter ``state``
        State with minimum description length.
47
48
49
50
51

    Notes
    -----

    This function is a convenience wrapper around
52
    :meth:`~graph_tool.inference.MultilevelMCMCState.multilevel_mcmc_sweep`.
53
54
55
56
57
58
59

    See [peixoto-efficient-2014]_ for details on the algorithm.

    This algorithm has a complexity of :math:`O(V \ln^2 V)`, where :math:`V` is
    the number of nodes in the network.

    Examples
Tiago Peixoto's avatar
Tiago Peixoto committed
60
    --------
61
62
63

    .. testsetup:: mdl

Tiago Peixoto's avatar
Tiago Peixoto committed
64
65
       gt.seed_rng(43)
       np.random.seed(43)
66
67
68
69
70
71

    .. doctest:: mdl

       >>> g = gt.collection.data["polbooks"]
       >>> state = gt.minimize_blockmodel_dl(g)
       >>> state.draw(pos=g.vp["pos"], vertex_shape=state.get_blocks(),
Tiago Peixoto's avatar
Tiago Peixoto committed
72
       ...            output="polbooks_blocks_mdl.svg")
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
       <...>

    .. figure:: polbooks_blocks_mdl.*
       :align: center

       Block partition of a political books network, which minimizes the
       description length of the network according to the degree-corrected
       stochastic blockmodel.


    .. testsetup:: mdl_overlap

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

    .. doctest:: mdl_overlap

       >>> g = gt.collection.data["polbooks"]
Tiago Peixoto's avatar
Tiago Peixoto committed
91
       >>> state = gt.minimize_blockmodel_dl(g, state=gt.OverlapBlockState)
Tiago Peixoto's avatar
Tiago Peixoto committed
92
       >>> state.draw(pos=g.vp["pos"], output="polbooks_overlap_blocks_mdl.svg")
93
94
95
96
97
98
99
100
101
       <...>

    .. figure:: polbooks_overlap_blocks_mdl.*
       :align: center

       Overlapping partition of a political books network, which minimizes the
       description length of the network according to the overlapping
       degree-corrected stochastic blockmodel.

Tiago Peixoto's avatar
Tiago Peixoto committed
102
103
104
105
106
    .. doctest:: mdl_pp

       >>> g = gt.collection.data["celegansneural"]
       >>> state = gt.minimize_blockmodel_dl(g, state=gt.PPBlockState)
       >>> state.draw(output="celegans_mdl_pp.pdf")
Tiago Peixoto's avatar
Tiago Peixoto committed
107
       <...>
Tiago Peixoto's avatar
Tiago Peixoto committed
108
109
110
111
112
113
114
115
116
117
118
119

    .. testcleanup:: mdl_pp

       conv_png("celegans_mdl_pp.pdf")

    .. figure:: celegans_mdl_pp.png
       :align: center
       :width: 60%

       Assortative partition of the *C. elegans* neural network, which minimizes
       the description length of the network according to the degree-corrected
       planted-partition blockmodel.
120
121
122
123
124
125
126
127

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

128
    state = state(g, **state_args)
129
130
131
132
133

    args = dict(niter=1, psingle=0, beta=numpy.inf)
    args.update(multilevel_mcmc_args)

    state.multilevel_mcmc_sweep(**args)
134
135
136

    return state

137
def minimize_nested_blockmodel_dl(g, init_bs=None,
138
                                  state=NestedBlockState, state_args={},
139
                                  multilevel_mcmc_args={}):
140
    r"""Fit the nested stochastic block model, by minimizing its description length
141
    using an agglomerative heuristic.
142
143
144
145
146

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        The graph.
147
148
149
    init_bs : iterable of iterable of ``int``s (optional, default: ``None``)
        Initial hierarchical partition.
    B_min : ``int`` (optional, default: ``1``)
150
        The minimum number of blocks.
151
    B_max : ``int`` (optional, default: ``numpy.iinfo(numpy.int64).max``)
152
        The maximum number of blocks.
153
    b_min : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
154
        The partition to be used with the minimum number of blocks.
155
    b_max : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
156
        The partition to be used with the maximum number of blocks.
157
    state : SBM state class (optional, default: :class:`~graph_tool.inference.NestedBlockState`)
158
        Type of model that will be used.
159
160
    state_args : ``dict`` (optional, default: ``{}``)
        Arguments to be passed to appropriate state constructor (e.g.
161
        :class:`~graph_tool.inference.BlockState`)
162
    multilevel_mcmc_args : ``dict`` (optional, default: ``{}``)
163
        Arguments to be passed to :meth:`~graph_tool.inference.MultilevelMCMCState.multilevel_mcmc_sweep`.
164
165
166

    Returns
    -------
167
168
    min_state : type given by parameter ``state``
        State with minimum description length.
169
170
171
172

    Notes
    -----
    This function is a convenience wrapper around
173
    :meth:`~graph_tool.inference.NestedBlockState.multilevel_mcmc_sweep`.
174
175
176

    See [peixoto-hierarchical-2014]_ for details on the algorithm.

Tiago Peixoto's avatar
Tiago Peixoto committed
177
178
    This algorithm has a complexity of :math:`O(E \ln^2 V)`, where :math:`E` and
    :math:`V` are the number of edges and nodes in the network, respectively.
179
180
181
182
183

    Examples
    --------
    .. testsetup:: nested_mdl

184
185
       gt.seed_rng(43)
       np.random.seed(43)
186
187
188
189

    .. doctest:: nested_mdl

       >>> g = gt.collection.data["power"]
Tiago Peixoto's avatar
Tiago Peixoto committed
190
       >>> state = gt.minimize_nested_blockmodel_dl(g)
191
192
193
194
195
       >>> state.draw(output="power_nested_mdl.pdf")
       (...)

    .. testcleanup:: nested_mdl

Tiago Peixoto's avatar
Tiago Peixoto committed
196
       conv_png("power_nested_mdl.pdf")
197

Tiago Peixoto's avatar
Tiago Peixoto committed
198
    .. figure:: power_nested_mdl.png
199
       :align: center
Tiago Peixoto's avatar
Tiago Peixoto committed
200
       :width: 60%
201
202
203
204
205
206
207
208
209

       Hierarchical Block partition of a power-grid network, which minimizes
       the description length of the network according to the nested
       (degree-corrected) stochastic blockmodel.


    .. doctest:: nested_mdl_overlap

       >>> g = gt.collection.data["celegansneural"]
Tiago Peixoto's avatar
Tiago Peixoto committed
210
       >>> state = gt.minimize_nested_blockmodel_dl(g, state_args=dict(overlap=True))
211
212
213
214
215
       >>> state.draw(output="celegans_nested_mdl_overlap.pdf")
       (...)

    .. testcleanup:: nested_mdl_overlap

Tiago Peixoto's avatar
Tiago Peixoto committed
216
       conv_png("celegans_nested_mdl_overlap.pdf")
217

Tiago Peixoto's avatar
Tiago Peixoto committed
218
    .. figure:: celegans_nested_mdl_overlap.png
219
       :align: center
Tiago Peixoto's avatar
Tiago Peixoto committed
220
       :width: 60%
221
222
223
224
225
226
227
228
229
230
231

       Overlapping block partition of the *C. elegans* neural network, which
       minimizes the description length of the network according to the nested
       overlapping degree-corrected stochastic blockmodel.

    References
    ----------
    .. [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`.
Tiago Peixoto's avatar
Tiago Peixoto committed
232

233
234
    """

235
    state = state(g, bs=init_bs, **state_args)
236
237
238
239
240
241
242
243
244
245
246
247
248
249

    args = dict(niter=1, psingle=0, beta=numpy.inf)
    args.update(multilevel_mcmc_args)

    l = 0
    while l >= 0:

        ret = state.multilevel_mcmc_sweep(ls=[l], **args)

        if args.get("verbose", False):
            print(l, ret, state)

        if abs(ret[0]) < 1e-8:
            l -= 1
250
        else:
251
            l = min(l + 1, len(state.levels) - 1)
252

253
    return state