Implement spectral functions in C++

configure.ac

@@ -403,18 +403,19 @@ AC_CONFIG_FILES([
 Makefile
 src/Makefile
 src/graph/Makefile
-src/graph/correlations/Makefile
-src/graph/generation/Makefile
-src/graph/stats/Makefile
 src/graph/centrality/Makefile
 src/graph/clustering/Makefile
 src/graph/community/Makefile
-src/graph/util/Makefile
-src/graph/topology/Makefile
+src/graph/correlations/Makefile
+src/graph/draw/Makefile
 src/graph/flow/Makefile
+src/graph/generation/Makefile
 src/graph/layout/Makefile
 src/graph/search/Makefile
-src/graph/draw/Makefile
+src/graph/spectral/Makefile
+src/graph/stats/Makefile
+src/graph/topology/Makefile
+src/graph/util/Makefile
 src/graph_tool/Makefile
 ])
 AC_OUTPUT

src/graph/Makefile.am

 ## Process this file with automake to produce Makefile.in
 
-SUBDIRS = centrality clustering community correlations draw flow generation layout search stats topology util
+SUBDIRS = centrality clustering community correlations draw flow generation layout search spectral stats topology util
 
 AM_CPPFLAGS =\
     -I$(srcdir)/.. \

src/graph/graph_selectors.hh

@@ -239,6 +239,21 @@ struct in_edge_iteratorS
     }
 };
 
+// out edges selector for completeness
+template <class Graph>
+struct out_edge_iteratorS
+{
+    typedef typename graph_traits<Graph>::out_edge_iterator type;
+
+    typedef typename boost::graph_traits<Graph>::vertex_descriptor
+        vertex_descriptor;
+    static std::pair<type,type> get_edges(vertex_descriptor v,
+                                          const Graph& g)
+    {
+        return out_edges(v, g);
+    }
+};
+
 // helper types for all_edges_iteratorS
 template <class Graph, class IsDirected>
 struct get_all_edges

src/graph/spectral/Makefile.am

+## Process this file with automake to produce Makefile.in
+
+AM_CPPFLAGS = $(MOD_CPPFLAGS) $(CAIROMM_CFLAGS)
+
+AM_CFLAGS = $(AM_CXXFLAGS)
+
+libgraph_tool_spectraldir = $(MOD_DIR)/spectral
+
+libgraph_tool_spectral_LTLIBRARIES = libgraph_tool_spectral.la
+
+libgraph_tool_spectral_la_includedir = $(MOD_DIR)/include
+
+libgraph_tool_spectral_la_LIBADD = $(MOD_LIBADD)
+
+libgraph_tool_spectral_la_LDFLAGS = $(MOD_LDFLAGS)
+
+libgraph_tool_spectral_la_SOURCES = \
+    graph_adjacency.cc \
+    graph_incidence.cc \
+    graph_laplacian.cc \
+    graph_norm_laplacian.cc \
+    graph_matrix.cc
+
+libgraph_tool_spectral_la_include_HEADERS = \
+    graph_adjacency.hh \
+    graph_incidence.hh \
+    graph_laplacian.hh

src/graph/spectral/graph_adjacency.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_adjacency.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void adjacency(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, bind<void>(get_adjacency(),
+                       _1, _2, _3, ref(data), ref(i), ref(j)),
+         vertex_scalar_properties(),
+         weight_props_t())(index, weight);
+
+}

src/graph/spectral/graph_adjacency.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_ADJACENCY_MATRIX_HH
+#define GRAPH_ADJACENCY_MATRIX_HH
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+
+struct get_adjacency
+{
+    template <class Graph, class Index, class Weight>
+    void operator()(Graph& g, Index index, Weight weight,
+                    multi_array_ref<double,1>& data,
+                    multi_array_ref<int,1>& i,
+                    multi_array_ref<int,1>& j) const
+    {
+        int pos = 0;
+        typename graph_traits<Graph>::edge_iterator e, e_end;
+        for(tie(e, e_end) = edges(g); e != e_end; ++e)
+        {
+            data[pos] = get(weight, *e);
+            i[pos] = get(index, target(*e, g));
+            j[pos] = get(index, source(*e, g));
+
+            ++pos;
+            if (!is_directed::apply<Graph>::type::value)
+            {
+                data[pos] = get(weight, *e);
+                i[pos] = get(index, source(*e, g));
+                j[pos] = get(index, target(*e, g));
+                ++pos;
+            }
+        }
+    }
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_ADJACENCY_MATRIX_HH

src/graph/spectral/graph_incidence.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_incidence.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void incidence(GraphInterface& g, boost::any vindex, boost::any eindex,
+               python::object odata, python::object oi,
+               python::object oj)
+{
+    if (!belongs<vertex_scalar_properties>()(vindex))
+        throw ValueException("index vertex property must have a scalar value type");
+    if (!belongs<edge_scalar_properties>()(eindex))
+        throw ValueException("index edge property must have a scalar value type");
+
+    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, bind<void>(get_incidence(),
+                       _1, _2, _3, ref(data), ref(i), ref(j)),
+         vertex_scalar_properties(),
+         edge_scalar_properties())(vindex, eindex);
+
+}

src/graph/spectral/graph_incidence.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_INCIDENCE_HH
+#define GRAPH_INCIDENCE_HH
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+
+struct get_incidence
+{
+    template <class Graph, class VIndex, class EIndex>
+    void operator()(Graph& g, VIndex vindex, EIndex eindex,
+                    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)
+        {
+
+            typename graph_traits<Graph>::out_edge_iterator e, e_end;
+            for(tie(e, e_end) = out_edges(*v, g); e != e_end; ++e)
+            {
+                if (is_directed::apply<Graph>::type::value)
+                    data[pos] = -1;
+                else
+                    data[pos] = 1;
+                i[pos] = get(vindex, *v);
+                j[pos] = get(eindex, *e);
+                ++pos;
+            }
+
+            typename in_edge_iteratorS<Graph>::type ie, ie_end;
+            for(tie(ie, ie_end) = in_edge_iteratorS<Graph>::get_edges(*v, g);
+                ie != ie_end; ++ie)
+            {
+                data[pos] = 1;
+                i[pos] = get(vindex, *v);
+                j[pos] = get(eindex, *ie);
+                ++pos;
+            }
+        }
+    }
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_INCIDENCE_HH

src/graph/spectral/graph_laplacian.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_laplacian.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void laplacian(GraphInterface& g, boost::any index, boost::any weight,
+               string sdeg,
+               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);
+
+    deg_t deg;
+    if (sdeg == "in")
+        deg = IN_DEG;
+    if (sdeg == "out")
+        deg = OUT_DEG;
+    if (sdeg == "total")
+        deg = TOTAL_DEG;
+
+    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, bind<void>(get_laplacian(),
+                       _1, _2, _3, deg, ref(data), ref(i), ref(j)),
+         vertex_scalar_properties(),
+         weight_props_t())(index, weight);
+
+}

src/graph/spectral/graph_laplacian.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_LAPLACIAN_HH
+#define GRAPH_LAPLACIAN_HH
+
+#include "graph.hh"
+#include "graph_filtering.hh"
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+
+enum deg_t
+{
+    IN_DEG,
+    OUT_DEG,
+    TOTAL_DEG
+};
+
+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, all_edges_iteratorS<Graph>)
+{
+    return total_degreeS()(v, g);
+}
+
+template <class Graph, class EdgeSelector>
+double
+sum_degree(Graph& g, typename graph_traits<Graph>::vertex_descriptor v,
+           ConstantPropertyMap<double, GraphInterface::edge_t> w, in_edge_iteratorS<Graph>)
+{
+    return in_degreeS()(v, g);
+}
+
+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_laplacian
+{
+    template <class Graph, class Index, class Weight>
+    void operator()(const Graph& g, Index index, Weight weight, deg_t deg,
+                    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>::edge_iterator e, e_end;
+        for(tie(e, e_end) = edges(g); e != e_end; ++e)
+        {
+            if (source(*e, g) == target(*e, g))
+                continue;
+
+            data[pos] = -get(weight, *e);
+            i[pos] = get(index, target(*e, g));
+            j[pos] = get(index, source(*e, g));
+
+            ++pos;
+            if (!is_directed::apply<Graph>::type::value)
+            {
+                data[pos] = -get(weight, *e);
+                i[pos] = get(index, source(*e, g));
+                j[pos] = get(index, target(*e, g));
+                ++pos;
+            }
+        }
+
+        typename graph_traits<Graph>::vertex_iterator v, v_end;
+        for(tie(v, v_end) = vertices(g); v != v_end; ++v)
+        {
+            double k = 0;
+            switch (deg)
+            {
+            case OUT_DEG:
+                k = sum_degree(g, *v, weight, out_edge_iteratorS<Graph>());
+                break;
+            case IN_DEG:
+                k = sum_degree(g, *v, weight, in_edge_iteratorS<Graph>());
+                break;
+            case TOTAL_DEG:
+                k = sum_degree(g, *v, weight, all_edges_iteratorS<Graph>());
+            }
+            data[pos] = k;
+            i[pos] = j[pos] = get(index, *v);
+            ++pos;
+        }
+
+    }
+};
+
+
+
+struct get_norm_laplacian
+{
+    template <class Graph, class Index, class Weight>
+    void operator()(const Graph& g, Index index, Weight weight, deg_t deg,
+                    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 ks = 0;
+            switch (deg)
+            {
+            case OUT_DEG:
+                ks = sum_degree(g, *v, weight, out_edge_iteratorS<Graph>());
+                break;
+            case IN_DEG:
+                ks = sum_degree(g, *v, weight, in_edge_iteratorS<Graph>());
+                break;
+            case TOTAL_DEG:
+                ks = sum_degree(g, *v, weight, all_edges_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)
+            {
+                if (source(*e, g) == target(*e, g))
+                    continue;
+                double kt = 0;
+
+                switch (deg)
+                {
+                case OUT_DEG:
+                    kt = sum_degree(g, target(*e, g), weight, out_edge_iteratorS<Graph>());
+                    break;
+                case IN_DEG:
+                    kt = sum_degree(g, target(*e, g), weight, in_edge_iteratorS<Graph>());
+                    break;
+                case TOTAL_DEG:
+                    kt = sum_degree(g, target(*e, g), weight, all_edges_iteratorS<Graph>());
+                }
+
+                if (ks * kt > 0)
+                    data[pos] = -get(weight, *e) / sqrt(ks * kt);
+                i[pos] = get(index, target(*e, g));
+                j[pos] = get(index, source(*e, g));
+
+                ++pos;
+            }
+
+            if (ks > 0)
+                data[pos] = 1;
+            i[pos] = j[pos] = get(index, *v);
+            ++pos;
+        }
+
+    }
+};
+
+
+} // namespace graph_tool
+
+#endif // GRAPH_LAPLACIAN_HH

src/graph/spectral/graph_matrix.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/>.
+
+#in