__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 : :mod:`~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, \\
            0 & \text{otherwise}
        \end{cases}

    In the case of weighted edges, the value 1 is replaced 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::

      gt.seed_rng(42)

    >>> g = gt.random_graph(100, lambda: (10, 10))
    >>> m = gt.adjacency(g)
    >>> print(m.todense())
    [[ 0.  0.  0. ...,  0.  1.  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.]
     [ 0.  0.  1. ...,  0.  0.  0.]]

    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)
    m = scipy.sparse.coo_matrix((data, (i,j)))
    m = m.tocsr()
    return m


@_limit_args({"deg": ["total", "in", "out"]})
def laplacian(g, deg="total", normalized=True, 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: True)
        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 : :mod:`~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::

      gt.seed_rng(42)

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

    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

    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)
    m = scipy.sparse.coo_matrix((data, (i,j)))
    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 : :mod:`~scipy.sparse.csr_matrix`
        The (sparse) adjacency 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. -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. -1. -1.]]

    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