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

"""
``graph_tool.spectral`` - Spectral properties
---------------------------------------------
24
25
26
27
28
29
30
31
32
33
34
35
36

Summary
+++++++

.. autosummary::
   :nosignatures:

   adjacency
   laplacian
   incidence

Contents
++++++++
Tiago Peixoto's avatar
Tiago Peixoto committed
37
38
"""

39
40
from __future__ import division, absolute_import, print_function

41
from .. import _degree, _prop, Graph, _limit_args
42
43
from .. stats import label_self_loops
import numpy
Tiago Peixoto's avatar
Tiago Peixoto committed
44
45
import scipy.sparse

46
47
from .. dl_import import dl_import
dl_import("from . import libgraph_tool_spectral")
Tiago Peixoto's avatar
Tiago Peixoto committed
48
49
50

__all__ = ["adjacency", "laplacian", "incidence"]

51

52
def adjacency(g, weight=None, index=None):
53
54
55
56
    r"""Return the adjacency matrix of the graph.

    Parameters
    ----------
57
    g : :class:`~graph_tool.Graph`
58
        Graph to be used.
59
    weight : :class:`~graph_tool.PropertyMap` (optional, default: True)
60
        Edge property map with the edge weights.
61
62
63
    index : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Vertex property map specifying the row/column indexes. If not provided, the
        internal vertex index is used.
64
65
66

    Returns
    -------
Tiago Peixoto's avatar
Tiago Peixoto committed
67
    a : :class:`~scipy.sparse.csr_matrix`
68
        The (sparse) adjacency matrix.
69
70
71
72
73
74
75
76
77
78

    Notes
    -----
    The adjacency matrix is defined as

    .. math::

        a_{i,j} =
        \begin{cases}
            1 & \text{if } v_i \text{ is adjacent to } v_j, \\
Tiago Peixoto's avatar
Tiago Peixoto committed
79
            2 & \text{if } i = j, \text{ the graph is undirected and there is a self-loop incident in } v_i, \\
80
81
82
            0 & \text{otherwise}
        \end{cases}

Tiago Peixoto's avatar
Tiago Peixoto committed
83
84
    In the case of weighted edges, the entry values are multiplied by the weight
    of the respective edge.
85

86
87
88
    In the case of networks with parallel edges, the entries in the matrix
    become simply the edge multiplicities.

89
90
    Examples
    --------
91
92
    .. testsetup::

Tiago Peixoto's avatar
Tiago Peixoto committed
93
94
       import scipy.linalg
       from pylab import *
95

Tiago Peixoto's avatar
Tiago Peixoto committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
    >>> g = gt.collection.data["polblogs"]
    >>> A = gt.adjacency(g)
    >>> ew, ev = scipy.linalg.eig(A.todense())

    >>> figure(figsize=(8, 2))
    <...>
    >>> scatter(real(ew), imag(ew), c=abs(ew))
    <...>
    >>> xlabel(r"$\operatorname{Re}(\lambda)$")
    <...>
    >>> ylabel(r"$\operatorname{Im}(\lambda)$")
    <...>
    >>> tight_layout()
    >>> savefig("adjacency-spectrum.pdf")

    .. testcode::
       :hide:

       savefig("adjacency-spectrum.png")

    .. figure:: adjacency-spectrum.*
        :align: center

        Adjacency matrix spectrum for the political blog network.
120
121
122

    References
    ----------
123
    .. [wikipedia-adjacency] http://en.wikipedia.org/wiki/Adjacency_matrix
124
125
    """

126
127
128
129
130
131
    if index is None:
        if g.get_vertex_filter()[0] != None:
            index = g.new_vertex_property("int64_t")
            index.fa = numpy.arange(g.num_vertices())
        else:
            index = g.vertex_index
Tiago Peixoto's avatar
Tiago Peixoto committed
132

133
134
135
136
    E = g.num_edges() if g.is_directed() else 2 * g.num_edges()
    data = numpy.zeros(E, dtype="double")
    i = numpy.zeros(E, dtype="int32")
    j = numpy.zeros(E, dtype="int32")
Tiago Peixoto's avatar
Tiago Peixoto committed
137

138
139
    libgraph_tool_spectral.adjacency(g._Graph__graph, _prop("v", g, index),
                                     _prop("e", g, weight), data, i, j)
Tiago Peixoto's avatar
Tiago Peixoto committed
140
141
    V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1))
    m = scipy.sparse.coo_matrix((data, (i,j)), shape=(V, V))
142
143
    m = m.tocsr()
    return m
Tiago Peixoto's avatar
Tiago Peixoto committed
144

Tiago Peixoto's avatar
Tiago Peixoto committed
145
146

@_limit_args({"deg": ["total", "in", "out"]})
Tiago Peixoto's avatar
Tiago Peixoto committed
147
def laplacian(g, deg="total", normalized=False, weight=None, index=None):
148
149
150
151
    r"""Return the Laplacian matrix of the graph.

    Parameters
    ----------
152
    g : :class:`~graph_tool.Graph`
153
154
155
        Graph to be used.
    deg : str (optional, default: "total")
        Degree to be used, in case of a directed graph.
Tiago Peixoto's avatar
Tiago Peixoto committed
156
    normalized : bool (optional, default: False)
157
        Whether to compute the normalized Laplacian.
158
    weight : :class:`~graph_tool.PropertyMap` (optional, default: True)
159
        Edge property map with the edge weights.
160
161
162
    index : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Vertex property map specifying the row/column indexes. If not provided, the
        internal vertex index is used.
163
164
165

    Returns
    -------
Tiago Peixoto's avatar
Tiago Peixoto committed
166
    l : :class:`~scipy.sparse.csr_matrix`
167
        The (sparse) Laplacian matrix.
168
169
170

    Notes
    -----
Tiago Peixoto's avatar
Tiago Peixoto committed
171
    The weighted Laplacian matrix is defined as
172
173
174

    .. math::

Tiago Peixoto's avatar
Tiago Peixoto committed
175
        \ell_{ij} =
176
177
        \begin{cases}
        \Gamma(v_i) & \text{if } i = j \\
Tiago Peixoto's avatar
Tiago Peixoto committed
178
        -w_{ij}     & \text{if } i \neq j \text{ and } v_i \text{ is adjacent to } v_j \\
179
180
181
        0           & \text{otherwise}.
        \end{cases}

182
183
    Where :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij}` is sum of the weights of the
    edges incident on vertex :math:`v_i`. The normalized version is
184
185
186

    .. math::

Tiago Peixoto's avatar
Tiago Peixoto committed
187
        \ell_{ij} =
188
189
        \begin{cases}
        1         & \text{ if } i = j \text{ and } \Gamma(v_i) \neq 0 \\
Tiago Peixoto's avatar
Tiago Peixoto committed
190
        -\frac{w_{ij}}{\sqrt{\Gamma(v_i)\Gamma(v_j)}} & \text{ if } i \neq j \text{ and } v_i \text{ is adjacent to } v_j \\
191
192
193
        0         & \text{otherwise}.
        \end{cases}

Tiago Peixoto's avatar
Tiago Peixoto committed
194
195
196
197
198
    In the case of unweighted edges, it is assumed :math:`w_{ij} = 1`.

    For directed graphs, it is assumed :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij} +
    \sum_j A_{ji}w_{ji}` if ``deg=="total"``, :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij}`
    if ``deg=="out"`` or :math:`\Gamma(v_i)=\sum_j A_{ji}w_{ji}` ``deg=="in"``.
199
200
201

    Examples
    --------
Tiago Peixoto's avatar
Tiago Peixoto committed
202

203
204
    .. testsetup::

Tiago Peixoto's avatar
Tiago Peixoto committed
205
206
       import scipy.linalg
       from pylab import *
207

Tiago Peixoto's avatar
Tiago Peixoto committed
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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
    >>> g = gt.collection.data["polblogs"]
    >>> L = gt.laplacian(g)
    >>> ew, ev = scipy.linalg.eig(L.todense())

    >>> figure(figsize=(8, 2))
    <...>
    >>> scatter(real(ew), imag(ew), c=abs(ew))
    <...>
    >>> xlabel(r"$\operatorname{Re}(\lambda)$")
    <...>
    >>> ylabel(r"$\operatorname{Im}(\lambda)$")
    <...>
    >>> tight_layout()
    >>> savefig("laplacian-spectrum.pdf")

    .. testcode::
       :hide:

       savefig("laplacian-spectrum.png")

    .. figure:: laplacian-spectrum.*
        :align: center

        Laplacian matrix spectrum for the political blog network.

    >>> L = gt.laplacian(g, normalized=True)
    >>> ew, ev = scipy.linalg.eig(L.todense())

    >>> figure(figsize=(8, 2))
    <...>
    >>> scatter(real(ew), imag(ew), c=abs(ew))
    <...>
    >>> xlabel(r"$\operatorname{Re}(\lambda)$")
    <...>
    >>> ylabel(r"$\operatorname{Im}(\lambda)$")
    <...>
    >>> tight_layout()
    >>> savefig("norm-laplacian-spectrum.pdf")

    .. testcode::
       :hide:

       savefig("norm-laplacian-spectrum.png")

    .. figure:: norm-laplacian-spectrum.*
        :align: center

        Normalized Laplacian matrix spectrum for the political blog network.
256
257
258

    References
    ----------
259
    .. [wikipedia-laplacian] http://en.wikipedia.org/wiki/Laplacian_matrix
260
261
    """

262
263
264
265
266
267
268
    if index is None:
        if g.get_vertex_filter()[0] != None:
            index = g.new_vertex_property("int64_t")
            index.fa = numpy.arange(g.num_vertices())
        else:
            index = g.vertex_index

Tiago Peixoto's avatar
Tiago Peixoto committed
269
    V = g.num_vertices()
270
271
272
273
274
275
276
277
278
279
280
281
    nself = label_self_loops(g, mark_only=True).a.sum()
    E = g.num_edges() - nself
    if not g.is_directed():
        E *= 2
    N = E + g.num_vertices()
    data = numpy.zeros(N, dtype="double")
    i = numpy.zeros(N, dtype="int32")
    j = numpy.zeros(N, dtype="int32")

    if normalized:
        libgraph_tool_spectral.norm_laplacian(g._Graph__graph, _prop("v", g, index),
                                              _prop("e", g, weight), deg, data, i, j)
Tiago Peixoto's avatar
Tiago Peixoto committed
282
    else:
283
284
        libgraph_tool_spectral.laplacian(g._Graph__graph, _prop("v", g, index),
                                         _prop("e", g, weight), deg, data, i, j)
Tiago Peixoto's avatar
Tiago Peixoto committed
285
286
    V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1))
    m = scipy.sparse.coo_matrix((data, (i, j)), shape=(V, V))
287
    m = m.tocsr()
Tiago Peixoto's avatar
Tiago Peixoto committed
288
289
    return m

Tiago Peixoto's avatar
Tiago Peixoto committed
290

291
def incidence(g, vindex=None, eindex=None):
292
293
294
295
    r"""Return the incidence matrix of the graph.

    Parameters
    ----------
296
    g : :class:`~graph_tool.Graph`
297
        Graph to be used.
298
299
300
301
302
303
    vindex : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Vertex property map specifying the row indexes. If not provided, the
        internal vertex index is used.
    eindex : :class:`~graph_tool.PropertyMap` (optional, default: None)
        Edge property map specifying the column indexes. If not provided, the
        internal edge index is used.
304
305
306

    Returns
    -------
Tiago Peixoto's avatar
Tiago Peixoto committed
307
308
    a : :class:`~scipy.sparse.csr_matrix`
        The (sparse) incidence matrix.
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

    Notes
    -----
    For undirected graphs, the incidence matrix is defined as

    .. math::

        b_{i,j} =
        \begin{cases}
            1 & \text{if vertex } v_i \text{and edge } e_j \text{ are incident}, \\
            0 & \text{otherwise}
        \end{cases}

    For directed graphs, the definition is

    .. math::

        b_{i,j} =
        \begin{cases}
            1 & \text{if edge } e_j \text{ enters vertex } v_i, \\
            -1 & \text{if edge } e_j \text{ leaves vertex } v_i, \\
            0 & \text{otherwise}
        \end{cases}

    Examples
    --------
335
336
337
338
    .. testsetup::

      gt.seed_rng(42)

339
340
    >>> g = gt.random_graph(100, lambda: (2,2))
    >>> m = gt.incidence(g)
341
    >>> print(m.todense())
Tiago Peixoto's avatar
Tiago Peixoto committed
342
    [[-1. -1.  0. ...,  0.  0.  0.]
Tiago Peixoto's avatar
Tiago Peixoto committed
343
     [ 0.  0.  0. ...,  0.  0.  0.]
344
     [ 0.  0.  0. ...,  0.  0.  0.]
Tiago Peixoto's avatar
Tiago Peixoto committed
345
     ..., 
Tiago Peixoto's avatar
Tiago Peixoto committed
346
     [ 0.  0. -1. ...,  0.  0.  0.]
347
     [ 0.  0.  0. ...,  0.  0.  0.]
Tiago Peixoto's avatar
Tiago Peixoto committed
348
     [ 0.  0.  0. ...,  1.  0.  0.]]
349
350
351

    References
    ----------
352
    .. [wikipedia-incidence] http://en.wikipedia.org/wiki/Incidence_matrix
353
354
    """

355
356
357
358
359
360
    if vindex is None:
        if g.get_edge_filter()[0] != None:
            vindex = g.new_vertex_property("int64_t")
            vindex.fa = numpy.arange(g.num_vertices())
        else:
            vindex = g.vertex_index
Tiago Peixoto's avatar
Tiago Peixoto committed
361

362
363
364
365
366
367
    if eindex is None:
        if g.get_edge_filter()[0] != None:
            eindex = g.new_edge_property("int64_t")
            eindex.fa = numpy.arange(g.num_edges())
        else:
            eindex = g.edge_index
Tiago Peixoto's avatar
Tiago Peixoto committed
368
369

    E = g.num_edges()
370
371
372
373
374
375
376
377
    data = numpy.zeros(2 * E, dtype="double")
    i = numpy.zeros(2 * E, dtype="int32")
    j = numpy.zeros(2 * E, dtype="int32")

    libgraph_tool_spectral.incidence(g._Graph__graph, _prop("v", g, vindex),
                                     _prop("e", g, eindex), data, i, j)
    m = scipy.sparse.coo_matrix((data, (i,j)))
    m = m.tocsr()
Tiago Peixoto's avatar
Tiago Peixoto committed
378
    return m