Include transition() and modularity_matrix() in spectral module

src/graph/spectral/Makefile.am

@@ -19,9 +19,11 @@ libgraph_tool_spectral_la_SOURCES = \
     graph_incidence.cc \
     graph_laplacian.cc \
     graph_norm_laplacian.cc \
-    graph_matrix.cc
+    graph_matrix.cc \
+    graph_transition.cc
 
 libgraph_tool_spectral_la_include_HEADERS = \
     graph_adjacency.hh \
     graph_incidence.hh \
-    graph_laplacian.hh
+    graph_laplacian.hh \
+    graph_transition.hh

src/graph/spectral/graph_matrix.cc

@@ -42,6 +42,10 @@ void incidence(GraphInterface& g, boost::any vindex, boost::any eindex,
                python::object odata, python::object oi,
                python::object oj);
 
+void transition(GraphInterface& g, boost::any index, boost::any weight,
+                python::object odata, python::object oi,
+                python::object oj);
+
 BOOST_PYTHON_MODULE(libgraph_tool_spectral)
 {
     using namespace boost::python;
@@ -49,4 +53,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_spectral)
     def("laplacian", &laplacian);
     def("norm_laplacian", &norm_laplacian);
     def("incidence", &incidence);
+    def("transition", &transition);
 }

src/graph/spectral/graph_transition.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// 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/>.
+
+#include <boost/python.hpp>
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+#include "numpy_bind.hh"
+
+#include "graph_selectors.hh"
+#include "graph_properties.hh"
+
+#include "graph_transition.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void transition(GraphInterface& g, boost::any index, boost::any weight,
+                python::object odata, python::object oi,
+                python::object oj)
+{
+    if (!belongs<vertex_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    typedef ConstantPropertyMap<double, GraphInterface::edge_t> weight_map_t;
+    typedef mpl::push_back<edge_scalar_properties, weight_map_t>::type
+        weight_props_t;
+
+    if (!weight.empty() && !belongs<edge_scalar_properties>()(weight))
+        throw ValueException("weight edge property must have a scalar value type");
+
+    if(weight.empty())
+        weight = weight_map_t(1.0);
+
+    multi_array_ref<double,1> data = get_array<double,1>(odata);
+    multi_array_ref<int32_t,1> i = get_array<int32_t,1>(oi);
+    multi_array_ref<int32_t,1> j = get_array<int32_t,1>(oj);
+    run_action<>()
+        (g, std::bind(get_transition(),
+                      placeholders::_1,  placeholders::_2,  placeholders::_3,
+                      std::ref(data), std::ref(i), std::ref(j)),
+         vertex_scalar_properties(),
+         weight_props_t())(index, weight);
+
+}

src/graph/spectral/graph_transition.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// 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/>.
+
+#ifndef GRAPH_TRANSITION_HH
+#define GRAPH_TRANSITION_HH
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+
+using namespace boost;
+
+template <class Graph, class Weight, class EdgeSelector>
+typename property_traits<Weight>::value_type
+sum_degree(Graph& g, typename graph_traits<Graph>::vertex_descriptor v,
+           Weight w, EdgeSelector)
+{
+    typename property_traits<Weight>::value_type sum = 0;
+    typename EdgeSelector::type e, e_end;
+    for(tie(e, e_end) = EdgeSelector::get_edges(v, g); e != e_end; ++e)
+        sum += get(w, *e);
+    return sum;
+}
+
+template <class Graph, class EdgeSelector>
+double
+sum_degree(Graph& g, typename graph_traits<Graph>::vertex_descriptor v,
+           ConstantPropertyMap<double, GraphInterface::edge_t> w, out_edge_iteratorS<Graph>)
+{
+    return out_degreeS()(v, g);
+}
+
+struct get_transition
+{
+    template <class Graph, class Index, class Weight>
+    void operator()(const Graph& g, Index index, Weight weight,
+                    multi_array_ref<double,1>& data,
+                    multi_array_ref<int32_t,1>& i,
+                    multi_array_ref<int32_t,1>& j) const
+    {
+        int pos = 0;
+        typename graph_traits<Graph>::vertex_iterator v, v_end;
+        for(tie(v, v_end) = vertices(g); v != v_end; ++v)
+        {
+            double k = sum_degree(g, *v, weight, out_edge_iteratorS<Graph>());
+            typename graph_traits<Graph>::out_edge_iterator e, e_end;
+            for(tie(e, e_end) = out_edges(*v, g); e != e_end; ++e)
+            {
+                data[pos] = 1. / k;
+                i[pos] = get(index, source(*e, g));
+                j[pos] = get(index, target(*e, g));
+                ++pos;
+            }
+        }
+    }
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_TRANSITION_HH

src/graph_tool/spectral/__init__.py

@@ -31,6 +31,8 @@ Summary
    adjacency
    laplacian
    incidence
+   transition
+   modularity_matrix
 
 Contents
 ++++++++
@@ -42,11 +44,12 @@ from .. import _degree, _prop, Graph, _limit_args
 from .. stats import label_self_loops
 import numpy
 import scipy.sparse
+import scipy.sparse.linalg
 
 from .. dl_import import dl_import
 dl_import("from . import libgraph_tool_spectral")
 
-__all__ = ["adjacency", "laplacian", "incidence"]
+__all__ = ["adjacency", "laplacian", "incidence", "transition", "modularity_matrix"]
 
 
 def adjacency(g, weight=None, index=None):
@@ -376,3 +379,197 @@ def incidence(g, vindex=None, eindex=None):
     m = scipy.sparse.coo_matrix((data, (i,j)))
     m = m.tocsr()
     return m
+
+def transition(g, weight=None, index=None):
+    r"""Return the transition 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
+    -------
+    T : :class:`~scipy.sparse.csr_matrix`
+        The (sparse) transition matrix.
+
+    Notes
+    -----
+    The transition matrix is defined as
+
+    .. math::
+
+        T_{ij} = \frac{A_{ij}}{k_i}
+
+    where :math:`k_i = \sum_j A_{ij}`, and :math:`A_{ij}` is the adjacency
+    matrix.
+
+    In the case of weighted edges, the values of the adjacency matrix are
+    multiplied by the edge weights.
+
+    Examples
+    --------
+    .. testsetup::
+
+       import scipy.linalg
+       from pylab import *
+
+    >>> g = gt.collection.data["polblogs"]
+    >>> T = gt.transition(g)
+    >>> ew, ev = scipy.linalg.eig(T.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("transition-spectrum.pdf")
+
+    .. testcode::
+       :hide:
+
+       savefig("transition-spectrum.png")
+
+    .. figure:: transition-spectrum.*
+        :align: center
+
+        Transition matrix spectrum for the political blog network.
+
+    References
+    ----------
+    .. [wikipedia-transition] https://en.wikipedia.org/wiki/Stochastic_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.transition(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
+
+
+
+def modularity_matrix(g, weight=None, index=None):
+    r"""Return the modularity 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
+    -------
+    B : :class:`~scipy.sparse.linalg.LinearOperator`
+        The (sparse) modularity matrix, represented as a
+        :class:`~scipy.sparse.linalg.LinearOperator`.
+
+    Notes
+    -----
+    The modularity matrix is defined as
+
+    .. math::
+
+        B_{ij} =  A_{ij} - \frac{k^+_i k^-_j}{2E}
+
+    where :math:`k^+_i = \sum_j A_{ij}`, :math:`k^-_i = \sum_j A_{ji}`,
+    :math:`2E=\sum_{ij}A_{ij}` and :math:`A_{ij}` is the adjacency matrix.
+
+    In the case of weighted edges, the values of the adjacency matrix are
+    multiplied by the edge weights.
+
+    Examples
+    --------
+
+    .. testsetup::
+
+       import scipy.linalg
+       from pylab import *
+
+    >>> g = gt.collection.data["polblogs"]
+    >>> B = gt.modularity_matrix(g)
+    >>> B = B * np.identity(B.shape[0])  # transform to a dense matrix
+    >>> ew, ev = scipy.linalg.eig(B)
+
+    >>> 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("modularity-spectrum.pdf")
+
+    .. testcode::
+       :hide:
+
+       savefig("modularity-spectrum.png")
+
+    .. figure:: modularity-spectrum.*
+        :align: center
+
+        Modularity matrix spectrum for the political blog network.
+
+    References
+    ----------
+    .. [newman-modularity]  M. E. J. Newman, M. Girvan, "Finding and evaluating
+       community structure in networks", Phys. Rev. E 69, 026113 (2004).
+       :doi:`10.1103/PhysRevE.69.026113`
+    """
+
+    A = adjacency(g, weight=weight, index=index)
+    if g.is_directed():
+        k_in = g.degree_property_map("in", weight=weight).fa
+    else:
+        k_in = g.degree_property_map("out", weight=weight).fa
+    k_out = g.degree_property_map("out", weight=weight).fa
+
+    N = A.shape[0]
+    E2 = float(k_out.sum())
+
+    def matvec(x):
+        M = x.shape[0]
+        if len(x.shape) > 1:
+            x = x.reshape(M)
+        nx = A * x - k_out * numpy.dot(k_in, x) / E2
+        return nx
+
+    def rmatvec(x):
+        M = x.shape[0]
+        if len(x.shape) > 1:
+            x = x.reshape(M)
+        nx = A.T * x - k_in * numpy.dot(k_out, x) / E2
+        return nx
+
+    B = scipy.sparse.linalg.LinearOperator((g.num_vertices(), g.num_vertices()),
+                                           matvec=matvec, rmatvec=rmatvec,
+                                           dtype="float")
+
+    return B