spectral: implement more efficient linear operator interface

src/graph/spectral/graph_adjacency.cc

@@ -61,3 +61,59 @@ void adjacency(GraphInterface& g, boost::any index, boost::any weight,
          },
          vertex_scalar_properties(), weight_props_t())(index, weight);
 }
+
+void adjacency_matvec(GraphInterface& g, boost::any index, boost::any weight,
+                      python::object ov, python::object oret)
+{
+    if (!belongs<vertex_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    typedef UnityPropertyMap<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();
+
+    multi_array_ref<double,1> v = get_array<double,1>(ov);
+    multi_array_ref<double,1> ret = get_array<double,1>(oret);
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& w)
+         {
+             return adj_matvec(graph, vi, w, v, ret);
+         },
+         vertex_scalar_properties(), weight_props_t())(index, weight);
+}
+
+void adjacency_matmat(GraphInterface& g, boost::any index, boost::any weight,
+                      python::object ov, python::object oret)
+{
+    if (!belongs<vertex_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    typedef UnityPropertyMap<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();
+
+    multi_array_ref<double,2> v = get_array<double,2>(ov);
+    multi_array_ref<double,2> ret = get_array<double,2>(oret);
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& w)
+         {
+             return adj_matmat(graph, vi, w, v, ret);
+         },
+         vertex_scalar_properties(), weight_props_t())(index, weight);
+}

src/graph/spectral/graph_adjacency.hh

@@ -53,6 +53,58 @@ struct get_adjacency
     }
 };
 
+template <class Graph, class Vindex, class Weight, class V>
+void adj_matvec(Graph& g, Vindex index, Weight w, V& x, V& ret)
+{
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             size_t i = get(index, v);
+             std::remove_reference_t<decltype(ret[i])> y = 0;
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                     y += get(w, e) * x[get(index, target(e, g))];
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                     y += x[get(index, u)];
+             }
+             ret[i] = y;
+         });
+}
+
+template <class Graph, class Vindex, class Weight, class M>
+void adj_matmat(Graph& g, Vindex index, Weight w, M& x, M& ret)
+{
+    size_t k = x.shape()[1];
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             size_t i = get(index, v);
+             auto y = ret[i];
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                 {
+                     auto w_e = get(w, e);
+                     for (size_t l = 0; l < k; ++l)
+                         y[l] += w_e * x[get(index, target(e, g))][l];
+                 }
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                     for (size_t l = 0; l < k; ++l)
+                         y[l] += x[get(index, u)][l];
+             }
+
+         });
+}
+
 } // namespace graph_tool
 
 #endif // GRAPH_ADJACENCY_MATRIX_HH

src/graph/spectral/graph_incidence.cc

@@ -53,3 +53,39 @@ void incidence(GraphInterface& g, boost::any vindex, boost::any eindex,
          },
          vertex_scalar_properties(), edge_scalar_properties())(vindex, eindex);
 }
+
+void incidence_matvec(GraphInterface& g, boost::any vindex, boost::any eindex,
+                      python::object ov, python::object oret, bool transpose)
+{
+    if (!belongs<vertex_scalar_properties>()(vindex))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    multi_array_ref<double,1> v = get_array<double,1>(ov);
+    multi_array_ref<double,1> ret = get_array<double,1>(oret);
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& ei)
+         {
+             return inc_matvec(graph, vi, ei, v, ret, transpose);
+         },
+         vertex_scalar_properties(), edge_scalar_properties())(vindex, eindex);
+}
+
+void incidence_matmat(GraphInterface& g, boost::any vindex, boost::any eindex,
+                      python::object ov, python::object oret, bool transpose)
+{
+    if (!belongs<vertex_scalar_properties>()(vindex))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    multi_array_ref<double,2> v = get_array<double,2>(ov);
+    multi_array_ref<double,2> ret = get_array<double,2>(oret);
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& ei)
+         {
+             return inc_matmat(graph, vi, ei, v, ret, transpose);
+         },
+         vertex_scalar_properties(), edge_scalar_properties())(vindex, eindex);
+}

src/graph/spectral/graph_incidence.hh

@@ -39,7 +39,7 @@ struct get_incidence
         {
             for (const auto& e : out_edges_range(v, g))
             {
-                if (graph_tool::is_directed(g))
+                if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
                     data[pos] = -1;
                 else
                     data[pos] = 1;
@@ -48,7 +48,7 @@ struct get_incidence
                 ++pos;
             }
 
-            if (graph_tool::is_directed(g))
+            if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
             {
                 for (const auto& e : in_edges_range(v, g))
                 {
@@ -62,6 +62,105 @@ struct get_incidence
     }
 };
 
+template <class Graph, class Vindex, class EIndex, class V>
+void inc_matvec(Graph& g, Vindex vindex, EIndex eindex, V& x, V& ret, bool transpose)
+{
+    if (!transpose)
+    {
+        parallel_vertex_loop
+            (g,
+             [&](auto v)
+             {
+                 auto& y = ret[get(vindex, v)];
+                 for (const auto& e : out_edges_range(v, g))
+                 {
+                    auto u = eindex[e];
+                    if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                        y -= x[u];
+                    else
+                        y += x[u];
+                 }
+
+                 if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                 {
+                    for (const auto& e : in_edges_range(v, g))
+                    {
+                        auto u = eindex[e];
+                        y += x[u];
+                    }
+                 }
+             });
+    }
+    else
+    {
+        parallel_edge_loop
+            (g,
+             [&](const auto& e)
+             {
+                 auto u = eindex[e];
+                 if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                     ret[u] = x[get(vindex, target(e, g))] - x[get(vindex, source(e, g))];
+                 else
+                     ret[u] = x[get(vindex, target(e, g))] + x[get(vindex, source(e, g))];
+             });
+    }
+}
+
+template <class Graph, class Vindex, class Eindex, class M>
+void inc_matmat(Graph& g, Vindex vindex, Eindex eindex, M& x, M& ret, bool transpose)
+{
+    size_t k = x.shape()[1];
+    if (!transpose)
+    {
+        parallel_vertex_loop
+            (g,
+             [&](auto v)
+             {
+                 auto y = ret[get(vindex, v)];
+                 for (const auto& e : out_edges_range(v, g))
+                 {
+                    auto u = eindex[e];
+                    for (size_t i = 0; i < k; ++i)
+                    {
+                        if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                            y[i] -= x[u][i];
+                        else
+                            y[i] += x[u][i];
+                    }
+                 }
+
+                 if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                 {
+                    for (const auto& e : in_edges_range(v, g))
+                    {
+                        auto u = eindex[e];
+                        for (size_t i = 0; i < k; ++i)
+                            y[i] += x[u][i];
+                    }
+                 }
+             });
+    }
+    else
+    {
+        parallel_edge_loop
+            (g,
+             [&](const auto& e)
+             {
+                 auto u = eindex[e];
+                 auto s = get(vindex, source(e, g));
+                 auto t = get(vindex, target(e, g));
+                 for (size_t i = 0; i < k; ++i)
+                 {
+                     if constexpr (graph_tool::is_directed_::apply<Graph>::type::value)
+                         ret[u][i] = x[t][i] - x[s][i];
+                     else
+                         ret[u][i] = x[t][i] + x[s][i];
+                 }
+             });
+    }
+}
+
+
 } // namespace graph_tool
 
 #endif // GRAPH_INCIDENCE_HH

src/graph/spectral/graph_laplacian.cc

@@ -70,3 +70,65 @@ void laplacian(GraphInterface& g, boost::any index, boost::any weight,
          },
          vertex_scalar_properties(), weight_props_t())(index, weight);
 }
+
+void laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight,
+                      boost::any deg, python::object ov, python::object oret)
+{
+    if (!belongs<vertex_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    typedef UnityPropertyMap<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();
+
+    multi_array_ref<double,1> v = get_array<double,1>(ov);
+    multi_array_ref<double,1> ret = get_array<double,1>(oret);
+
+    typedef typename vprop_map_t<double>::type deg_t;
+    deg_t::unchecked_t d = any_cast<deg_t>(deg).get_unchecked();
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& w)
+         {
+             return lap_matvec(graph, vi, w, d, v, ret);
+         },
+         vertex_scalar_properties(), weight_props_t())(index, weight);
+}
+
+void laplacian_matmat(GraphInterface& g, boost::any index, boost::any weight,
+                      boost::any deg, python::object ov, python::object oret)
+{
+    if (!belongs<vertex_scalar_properties>()(index))
+        throw ValueException("index vertex property must have a scalar value type");
+
+    typedef UnityPropertyMap<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();
+
+    typedef typename vprop_map_t<double>::type deg_t;
+    deg_t::unchecked_t d = any_cast<deg_t>(deg).get_unchecked();
+
+    multi_array_ref<double,2> v = get_array<double,2>(ov);
+    multi_array_ref<double,2> ret = get_array<double,2>(oret);
+
+    run_action<>()
+        (g,
+         [&](auto&& graph, auto&& vi, auto&& w)
+         {
+             return lap_matmat(graph, vi, w, d, v, ret);
+         },
+         vertex_scalar_properties(), weight_props_t())(index, weight);
+}

src/graph/spectral/graph_laplacian.hh

@@ -119,6 +119,76 @@ struct get_laplacian
     }
 };
 
+template <class Graph, class Vindex, class Weight, class Deg, class V>
+void lap_matvec(Graph& g, Vindex index, Weight w, Deg d, V& x, V& ret)
+{
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             std::remove_reference_t<decltype(ret[v])> y = 0;
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                 {
+                     auto u = target(e, g);
+                     if (u == v)
+                         continue;
+                     auto w_e = get(w, e);
+                     y += w_e * x[get(index, u)];
+                 }
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                 {
+                     if (u == v)
+                         continue;
+                     y += x[get(index, u)];
+                 }
+             }
+             ret[get(index, v)] = d[v] * x[get(index, v)] - y;
+         });
+}
+
+template <class Graph, class Vindex, class Weight, class Deg, class Mat>
+void lap_matmat(Graph& g, Vindex index, Weight w, Deg d, Mat& x, Mat& ret)
+{
+    size_t M = x.shape()[1];
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             auto vi = get(index, v);
+             auto y = ret[vi];
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                 {
+                     auto u = target(e, g);
+                     if (u == v)
+                         continue;
+                     auto w_e = get(w, e);
+                     auto ui = get(index, u);
+                     for (size_t i = 0; i < M; ++i)
+                         y[i] += w_e * x[ui][i];
+                 }
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                 {
+                     if (u == v)
+                         continue;
+                     auto ui = get(index, u);
+                     for (size_t i = 0; i < M; ++i)
+                         y[i] += x[ui][i];
+                 }
+             }
+             for (size_t i = 0; i < M; ++i)
+                 ret[vi][i] = d[v] * x[vi][i] - y[i];
+         });
+}
 
 
 struct get_norm_laplacian
@@ -130,41 +200,34 @@ struct get_norm_laplacian
                     multi_array_ref<int32_t,1>& j) const
     {
         int pos = 0;
+        std::vector<double> degs(num_vertices(g));
         for (auto v : vertices_range(g))
         {
-            double ks = 0;
+            double k = 0;
             switch (deg)
             {
             case OUT_DEG:
-                ks = sum_degree(g, v, weight, out_edge_iteratorS<Graph>());
+                k = sum_degree(g, v, weight, out_edge_iteratorS<Graph>());
                 break;
             case IN_DEG:
-                ks = sum_degree(g, v, weight, in_edge_iteratorS<Graph>());
+                k = sum_degree(g, v, weight, in_edge_iteratorS<Graph>());
                 break;
             case TOTAL_DEG:
-                ks = sum_degree(g, v, weight, all_edges_iteratorS<Graph>());
+                k = sum_degree(g, v, weight, all_edges_iteratorS<Graph>());
             }
+            degs[v] = sqrt(k);
+        }
 
+        for (auto v : vertices_range(g))
+        {
+            double ks = degs[v];
             for(const auto& e : out_edges_range(v, g))
             {
                 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>());
-                }
-
+                double kt = degs[target(e, g)];
                 if (ks * kt > 0)
-                    data[pos] = -get(weight, e) / sqrt(ks * kt);
+                    data[pos] = -get(weight, e) / (ks * kt);
                 i[pos] = get(index, target(e, g));
                 j[pos] = get(index, source(e, g));
 
@@ -180,6 +243,83 @@ struct get_norm_laplacian
     }
 };
 
+template <class Graph, class Vindex, class Weight, class Deg, class V>
+void nlap_matvec(Graph& g, Vindex index, Weight w, Deg id, V& x, V& ret)
+{
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             auto vi = get(index, v);
+             std::remove_reference_t<decltype(ret[v])> y = 0;
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                 {
+                     auto u = target(e, g);
+                     if (u == v)
+                         continue;
+                     auto w_e = get(w, e);
+                     y += w_e * x[get(index, u)] * id[u];
+                 }
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                 {
+                     if (u == v)
+                         continue;
+                     y += x[get(index, u)] * id[u];
+                 }
+             }
+             if (id[v] > 0)
+                 ret[vi] = x[vi] - y * id[v];
+         });
+}
+
+template <class Graph, class Vindex, class Weight, class Deg, class Mat>
+void nlap_matmat(Graph& g, Vindex index, Weight w, Deg id, Mat& x, Mat& ret)
+{
+    size_t M = x.shape()[1];
+    parallel_vertex_loop
+        (g,
+         [&](auto v)
+         {
+             auto vi = get(index, v);
+             auto y = ret[vi];
+             if constexpr (!std::is_same_v<Weight, UnityPropertyMap<double, GraphInterface::edge_t>>)
+             {
+                 for (auto e : in_or_out_edges_range(v, g))
+                 {
+                     auto u = target(e, g);
+                     if (u == v)
+                         continue;
+                     auto w_e = get(w, e);
+                     auto ui = get(index, u);
+                     for (size_t i = 0; i < M; ++i)
+                         y[i] += w_e * x[ui][i] * id[u];
+                 }
+             }
+             else
+             {
+                 for (auto u : in_or_out_neighbors_range(v, g))
+                 {
+                     if (u == v)
+                         continue;
+                     auto ui = get(index, u);
+                     for (size_t i = 0; i < M; ++i)
+                         y[i] += x[ui][i] * id[u];
+                 }
+             }
+
+             if (id[v] > 0)
+             {
+                 for (size_t i = 0; i < M; ++i)
+                     y[i] = x[vi][i] - y[i] * id[v];
+             }
+         });
+}
+
 
 } // namespace graph_tool
 

src/graph/spectral/graph_matrix.cc

@@ -26,41 +26,101 @@ void adjacency(GraphInterface& g, boost::any index, boost::any weight,
                python::object odata, python::object oi,
                python::object oj);
 
+void adjacency_matvec(GraphInterface& g, boost::any index, boost::any weight,
+                      python::object ov, python::object oret);
+
+void adjacency_matmat(GraphInterface& g, boost::any index, boost::any weight,
+                      python::object ov, python::object oret);
 
 void laplacian(GraphInterface& g, boost::any index, boost::any weight,
                string sdeg,
                python::object odata, python::object oi,
                python::object oj);
 
+void laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight,
+                      boost::any deg, python::object ov, python::object oret);
+
+void laplacian_matmat(GraphInterface