spectral: Add hashimoto() implementation

src/graph/spectral/Makefile.am

@@ -20,10 +20,12 @@ libgraph_tool_spectral_la_SOURCES = \
     graph_laplacian.cc \
     graph_norm_laplacian.cc \
     graph_matrix.cc \
-    graph_transition.cc
+    graph_transition.cc \
+    graph_nonbacktracking.cc
 
 libgraph_tool_spectral_la_include_HEADERS = \
     graph_adjacency.hh \
     graph_incidence.hh \
     graph_laplacian.hh \
-    graph_transition.hh
+    graph_transition.hh \
+    graph_nonbacktracking.hh

src/graph/spectral/graph_matrix.cc

@@ -46,6 +46,9 @@ void transition(GraphInterface& g, boost::any index, boost::any weight,
                 python::object odata, python::object oi,
                 python::object oj);
 
+void nonbacktracking(GraphInterface& gi, boost::any index,
+                     std::vector<int64_t>& i, std::vector<int64_t>& j);
+
 BOOST_PYTHON_MODULE(libgraph_tool_spectral)
 {
     using namespace boost::python;
@@ -55,4 +58,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_spectral)
     def("norm_laplacian", &norm_laplacian);
     def("incidence", &incidence);
     def("transition", &transition);
+    def("nonbacktracking", &nonbacktracking);
 }

src/graph/spectral/graph_nonbacktracking.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2018 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_nonbacktracking.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void nonbacktracking(GraphInterface& gi, boost::any index,
+                     std::vector<int64_t>& i, std::vector<int64_t>& j)
+{
+    if (!belongs<edge_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    run_action<>()
+        (gi, [&](auto& g, auto idx){ get_nonbacktracking(g, idx, i, j);},
+         edge_scalar_properties())(index);
+
+}

src/graph/spectral/graph_nonbacktracking.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2018 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_NONBACKTRACKING_MATRIX_HH
+#define GRAPH_NONBACKTRACKING_MATRIX_HH
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+using namespace boost;
+
+template <class Graph, class Index>
+void get_nonbacktracking(Graph& g, Index index,
+                         std::vector<int64_t>& i,
+                         std::vector<int64_t>& j)
+{
+    for (auto u : vertices_range(g))
+    {
+        for (auto e1 : out_edges_range(u, g))
+        {
+            auto v = target(e1, g);
+            int64_t idx1 = index[e1];
+            if (!graph_tool::is_directed(g))
+                idx1 = (idx1 << 1) + (u > v);
+            for (auto e2 : out_edges_range(v, g))
+            {
+                auto w = target(e2, g);
+                if (w == u)
+                    continue;
+                int64_t idx2 = index[e2];
+                if (!graph_tool::is_directed(g))
+                    idx2 = (idx2 << 1) + (v > w);
+                i.push_back(idx1);
+                j.push_back(idx2);
+            }
+        }
+    }
+}
+
+} // namespace graph_tool
+
+#endif // GRAPH_NONBACKTRACKING_MATRIX_HH

src/graph_tool/spectral/__init__.py

@@ -33,6 +33,7 @@ Summary
    incidence
    transition
    modularity_matrix
+   hashimoto
 
 Contents
 ++++++++
@@ -40,7 +41,7 @@ Contents
 
 from __future__ import division, absolute_import, print_function
 
-from .. import _degree, _prop, Graph, GraphView, _limit_args
+from .. import _degree, _prop, Graph, GraphView, _limit_args, Vector_int64_t
 from .. stats import label_self_loops
 import numpy
 import scipy.sparse
@@ -49,7 +50,8 @@ import scipy.sparse.linalg
 from .. dl_import import dl_import
 dl_import("from . import libgraph_tool_spectral")
 
-__all__ = ["adjacency", "laplacian", "incidence", "transition", "modularity_matrix"]
+__all__ = ["adjacency", "laplacian", "incidence", "transition",
+           "modularity_matrix", "hashimoto"]
 
 
 def adjacency(g, weight=None, index=None):
@@ -616,3 +618,97 @@ def modularity_matrix(g, weight=None, index=None):
                                            dtype="float")
 
     return B
+
+def hashimoto(g, weight=None, index=None):
+    r"""Return the Hashimoto (or non-backtracking) matrix of a graph.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    index : :class:`~graph_tool.PropertyMap` (optional, default: None)
+        Edge property map specifying the row/column indexes. If not provided, the
+        internal edge index is used.
+
+    Returns
+    -------
+    H : :class:`~scipy.sparse.csr_matrix`
+        The (sparse) Hashimoto matrix.
+
+    Notes
+    -----
+    The Hashimoto (a.k.a. non-backtracking) matrix is defined as
+
+    .. math::
+
+        h_{k\to l,i\to j} =
+        \begin{cases}
+            1 & \text{if } (k,l) \in E, (i,j) \in E, l=i, k\ne j,\\
+            0 & \text{otherwise},
+        \end{cases}
+
+    where :math:`E` is the edge set. It is therefore a :math:`2|E|\times 2|E|`
+    asymmetric square matrix (or :math:`|E|\times |E|` for directed graphs),
+    indexed over edge directions.
+
+    Examples
+    --------
+    .. testsetup::
+
+       import scipy.linalg
+       from pylab import *
+
+    >>> g = gt.collection.data["football"]
+    >>> H = gt.hashimoto(g)
+    >>> ew, ev = scipy.linalg.eig(H.todense())
+    >>> figure(figsize=(8, 4))
+    <...>
+    >>> scatter(real(ew), imag(ew), c=sqrt(abs(ew)), linewidths=0, alpha=0.6)
+    <...>
+    >>> xlabel(r"$\operatorname{Re}(\lambda)$")
+    Text(...)
+    >>> ylabel(r"$\operatorname{Im}(\lambda)$")
+    Text(...)
+    >>> tight_layout()
+    >>> savefig("hashimoto-spectrum.pdf")
+
+    .. testcode::
+       :hide:
+
+       savefig("hashimoto-spectrum.png")
+
+    .. figure:: hashimoto-spectrum.*
+        :align: center
+
+        Hashimoto matrix spectrum for the network of American football teams.
+
+    References
+    ----------
+    .. [hashimoto] Hashimoto, Ki-ichiro. "Zeta functions of finite graphs and
+       representations of p-adic groups." Automorphic forms and geometry of
+       arithmetic varieties. 1989. 211-280. :DOI:`10.1016/B978-0-12-330580-0.50015-X`
+
+    """
+
+    if index is None:
+        if g.get_edge_filter()[0] is not None:
+            index = g.new_edge_property("int64_t")
+            index.fa = numpy.arange(g.num_edges())
+            E = index.fa.max() + 1
+        else:
+            index = g.edge_index
+            E = g.edge_index_range
+
+    if not g.is_directed():
+        E *= 2
+
+    i = Vector_int64_t()
+    j = Vector_int64_t()
+
+    libgraph_tool_spectral.nonbacktracking(g._Graph__graph, _prop("e", g, index),
+                                           i, j)
+
+    data = numpy.ones(i.a.shape)
+    m = scipy.sparse.coo_matrix((data, (i.a,j.a)), shape=(E, E))
+    m = m.tocsr()
+    return m