__init__.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-
#
# graph_tool -- a general graph manipulation python module
#
# Copyright (C) 2006-2013 Tiago de Paula Peixoto <tiago@skewed.de>
#
# 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
---------------------------------------------

Summary
+++++++

.. autosummary::
   :nosignatures:

   adjacency
   laplacian
   incidence

Contents
++++++++
"""

from __future__ import division, absolute_import, print_function

from .. import _degree, _prop, Graph, _limit_args
from .. stats import label_self_loops
import numpy
import scipy.sparse

from .. dl_import import dl_import
dl_import("from . import libgraph_tool_spectral")

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


def adjacency(g, weight=None, index=None):
    r"""Return the adjacency matrix of the graph.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    weight : :class:`~graph_tool.PropertyMap` (optional, default: True)
        Edge property map with the edge weights.
    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.

    Returns
    -------
    a : :class:`~scipy.sparse.csr_matrix`
        The (sparse) adjacency matrix.

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

    .. math::

        a_{i,j} =
        \begin{cases}
            1 & \text{if } v_i \text{ is adjacent to } v_j, \\
            2 & \text{if } i = j, \text{ the graph is undirected and there is a self-loop incident in } v_i, \\
            0 & \text{otherwise}
        \end{cases}

    In the case of weighted edges, the entry values are multiplied by the weight
    of the respective edge.

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

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

       import scipy.linalg
       from pylab import *

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

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

    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

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

    libgraph_tool_spectral.adjacency(g._Graph__graph, _prop("v", g, index),
                                     _prop("e", g, weight), data, i, j)
    V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1))
    m = scipy.sparse.coo_matrix((data, (i,j)), shape=(V, V))
    m = m.tocsr()
    return m


@_limit_args({"deg": ["total", "in", "out"]})
def laplacian(g, deg="total", normalized=False, weight=None, index=None):
    r"""Return the Laplacian matrix of the graph.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    deg : str (optional, default: "total")
        Degree to be used, in case of a directed graph.
    normalized : bool (optional, default: False)
        Whether to compute the normalized Laplacian.
    weight : :class:`~graph_tool.PropertyMap` (optional, default: True)
        Edge property map with the edge weights.
    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.

    Returns
    -------
    l : :class:`~scipy.sparse.csr_matrix`
        The (sparse) Laplacian matrix.

    Notes
    -----
    The weighted Laplacian matrix is defined as

    .. math::

        \ell_{ij} =
        \begin{cases}
        \Gamma(v_i) & \text{if } i = j \\
        -w_{ij}     & \text{if } i \neq j \text{ and } v_i \text{ is adjacent to } v_j \\
        0           & \text{otherwise}.
        \end{cases}

    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

    .. math::

        \ell_{ij} =
        \begin{cases}
        1         & \text{ if } i = j \text{ and } \Gamma(v_i) \neq 0 \\
        -\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 \\
        0         & \text{otherwise}.
        \end{cases}

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

    Examples
    --------

    .. testsetup::

       import scipy.linalg
       from pylab import *

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

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

    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

    V = g.num_vertices()
    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)
    else:
        libgraph_tool_spectral.laplacian(g._Graph__graph, _prop("v", g, index),
                                         _prop("e", g, weight), deg, data, i, j)
    V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1))
    m = scipy.sparse.coo_matrix((data, (i, j)), shape=(V, V))
    m = m.tocsr()
    return m


def incidence(g, vindex=None, eindex=None):
    r"""Return the incidence matrix of the graph.

    Parameters
    ----------
    g : :class:`~graph_tool.Graph`
        Graph to be used.
    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.

    Returns
    -------
    a : :class:`~scipy.sparse.csr_matrix`
        The (sparse) incidence matrix.

    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
    --------
    .. testsetup::

      gt.seed_rng(42)

    >>> g = gt.random_graph(100, lambda: (2,2))
    >>> m = gt.incidence(g)
    >>> print(m.todense())
    [[-1. -1.  0. ...,  0.  0.  0.]
     [ 0.  0.  0. ...,  0.  0.  0.]
     [ 0.  0.  0. ...,  0.  0.  0.]
     ..., 
     [ 0.  0. -1. ...,  0.  0.  0.]
     [ 0.  0.  0. ...,  0.  0.  0.]
     [ 0.  0.  0. ...,  1.  0.  0.]]

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

    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

    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

    E = g.num_edges()
    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()
    return m