diff --git a/src/graph/spectral/graph_adjacency.cc b/src/graph/spectral/graph_adjacency.cc index 7a6896409e4a671e4aa5630e22f0c9db939f1474..d7f55ebd75a3a78eb184b49dbe316c947b9670ec 100644 --- a/src/graph/spectral/graph_adjacency.cc +++ b/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()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(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()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(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); +} diff --git a/src/graph/spectral/graph_adjacency.hh b/src/graph/spectral/graph_adjacency.hh index 469797c66e53004303f10b50c982a6f883b04710..b5552e109b2132a49fb2db243c388baad4552850 100644 --- a/src/graph/spectral/graph_adjacency.hh +++ b/src/graph/spectral/graph_adjacency.hh @@ -53,6 +53,58 @@ struct get_adjacency } }; +template +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 y = 0; + if constexpr (!std::is_same_v>) + { + 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 +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>) + { + 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 diff --git a/src/graph/spectral/graph_incidence.cc b/src/graph/spectral/graph_incidence.cc index 41930768705b84451b1845a430a568718cc04bf4..39ce49866a343420f5c4f39532a7c64ce7dc61d7 100644 --- a/src/graph/spectral/graph_incidence.cc +++ b/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()(vindex)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(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()(vindex)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(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); +} diff --git a/src/graph/spectral/graph_incidence.hh b/src/graph/spectral/graph_incidence.hh index b32d56af4183ac633ba178a227c41638db25d5a4..acb2c101b1a0c8e417fee5c54a7f91ea9dd92a1a 100644 --- a/src/graph/spectral/graph_incidence.hh +++ b/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::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::type::value) { for (const auto& e : in_edges_range(v, g)) { @@ -62,6 +62,105 @@ struct get_incidence } }; +template +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::type::value) + y -= x[u]; + else + y += x[u]; + } + + if constexpr (graph_tool::is_directed_::apply::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::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 +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::type::value) + y[i] -= x[u][i]; + else + y[i] += x[u][i]; + } + } + + if constexpr (graph_tool::is_directed_::apply::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::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 diff --git a/src/graph/spectral/graph_laplacian.cc b/src/graph/spectral/graph_laplacian.cc index 39d070ff402ebb6b56598912d267631b869bad32..e41db345fb668797de822ef0d74f33457d23bcec 100644 --- a/src/graph/spectral/graph_laplacian.cc +++ b/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()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(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()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(deg).get_unchecked(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(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); +} diff --git a/src/graph/spectral/graph_laplacian.hh b/src/graph/spectral/graph_laplacian.hh index d96ec94e3ea158e082d964608d89bdd3a85f9d3d..e08aa4845b5966c94b8b5438679f54b2e6ca86ba 100644 --- a/src/graph/spectral/graph_laplacian.hh +++ b/src/graph/spectral/graph_laplacian.hh @@ -119,6 +119,76 @@ struct get_laplacian } }; +template +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 y = 0; + if constexpr (!std::is_same_v>) + { + 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 +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>) + { + 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& j) const { int pos = 0; + std::vector 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()); + k = sum_degree(g, v, weight, out_edge_iteratorS()); break; case IN_DEG: - ks = sum_degree(g, v, weight, in_edge_iteratorS()); + k = sum_degree(g, v, weight, in_edge_iteratorS()); break; case TOTAL_DEG: - ks = sum_degree(g, v, weight, all_edges_iteratorS()); + k = sum_degree(g, v, weight, all_edges_iteratorS()); } + 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()); - break; - case IN_DEG: - kt = sum_degree(g, target(e, g), weight, in_edge_iteratorS()); - break; - case TOTAL_DEG: - kt = sum_degree(g, target(e, g), weight, all_edges_iteratorS()); - } - + 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 +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 y = 0; + if constexpr (!std::is_same_v>) + { + 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 +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>) + { + 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 diff --git a/src/graph/spectral/graph_matrix.cc b/src/graph/spectral/graph_matrix.cc index 8ba49da94fd55360eed1e1433945df12e485ff95..171da8c10f4ba146ab12d35912ffe7b6d8d3d78e 100644 --- a/src/graph/spectral/graph_matrix.cc +++ b/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& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret); void norm_laplacian(GraphInterface& g, boost::any index, boost::any weight, string sdeg, python::object odata, python::object oi, python::object oj); +void norm_laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret); + +void norm_laplacian_matmat(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret); + void incidence(GraphInterface& g, boost::any vindex, boost::any eindex, python::object odata, python::object oi, python::object oj); +void incidence_matvec(GraphInterface& g, boost::any vindex, boost::any eindex, + python::object ov, python::object oret, bool transpose); + +void incidence_matmat(GraphInterface& g, boost::any vindex, boost::any eindex, + python::object ov, python::object oret, bool tranpose); + void transition(GraphInterface& g, boost::any index, boost::any weight, python::object odata, python::object oi, python::object oj); +void transition_matvec(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret, + bool transpose); + +void transition_matmat(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret, + bool transpose); + void nonbacktracking(GraphInterface& gi, boost::any index, std::vector& i, std::vector& j); +void nonbacktracking_matvec(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose); + +void nonbacktracking_matmat(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose); + void compact_nonbacktracking(GraphInterface& gi, std::vector& i, std::vector& j, std::vector& x); +void compact_nonbacktracking_matvec(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose); + +void compact_nonbacktracking_matmat(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose); + BOOST_PYTHON_MODULE(libgraph_tool_spectral) { using namespace boost::python; docstring_options dopt(true, false); def("adjacency", &adjacency); + def("adjacency_matvec", &adjacency_matvec); + def("adjacency_matmat", &adjacency_matmat); def("laplacian", &laplacian); + def("laplacian_matvec", &laplacian_matvec); + def("laplacian_matmat", &laplacian_matmat); def("norm_laplacian", &norm_laplacian); + def("norm_laplacian_matvec", &norm_laplacian_matvec); + def("norm_laplacian_matmat", &norm_laplacian_matmat); def("incidence", &incidence); + def("incidence_matvec", &incidence_matvec); + def("incidence_matmat", &incidence_matmat); def("transition", &transition); + def("transition_matvec", &transition_matvec); + def("transition_matmat", &transition_matmat); def("nonbacktracking", &nonbacktracking); + def("nonbacktracking_matvec", &nonbacktracking_matvec); + def("nonbacktracking_matmat", &nonbacktracking_matmat); def("compact_nonbacktracking", &compact_nonbacktracking); + def("compact_nonbacktracking_matvec", &compact_nonbacktracking_matvec); + def("compact_nonbacktracking_matmat", &compact_nonbacktracking_matmat); } diff --git a/src/graph/spectral/graph_nonbacktracking.cc b/src/graph/spectral/graph_nonbacktracking.cc index c06bb469f96ed74819857bbb3afda5bd4c42ae4d..3f7eae93f1d0747aa024783bff48151ee5e241b2 100644 --- a/src/graph/spectral/graph_nonbacktracking.cc +++ b/src/graph/spectral/graph_nonbacktracking.cc @@ -42,6 +42,50 @@ void nonbacktracking(GraphInterface& gi, boost::any index, } +void nonbacktracking_matvec(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& ei) + { + if (!transpose) + nbt_matvec(graph, ei, v, ret); + else + nbt_matvec(graph, ei, v, ret); + }, + edge_scalar_properties())(index); +} + +void nonbacktracking_matmat(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& ei) + { + if (!transpose) + nbt_matmat(graph, ei, v, ret); + else + nbt_matmat(graph, ei, v, ret); + }, + edge_scalar_properties())(index); +} + void compact_nonbacktracking(GraphInterface& gi, std::vector& i, std::vector& j, std::vector& x) { @@ -49,3 +93,47 @@ void compact_nonbacktracking(GraphInterface& gi, std::vector& i, (gi, [&](auto& g){ get_compact_nonbacktracking(g, i, j, x);})(); } + +void compact_nonbacktracking_matvec(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& ei) + { + if (!transpose) + cnbt_matvec(graph, ei, v, ret); + else + cnbt_matvec(graph, ei, v, ret); + }, + vertex_scalar_properties())(index); +} + +void compact_nonbacktracking_matmat(GraphInterface& g, boost::any index, + python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& ei) + { + if (!transpose) + cnbt_matmat(graph, ei, v, ret); + else + cnbt_matmat(graph, ei, v, ret); + }, + vertex_scalar_properties())(index); +} diff --git a/src/graph/spectral/graph_nonbacktracking.hh b/src/graph/spectral/graph_nonbacktracking.hh index b7dd8a1d78ea03cc10e7c0820456bff891f63847..a263d63d6ad4f26ed84f5be2dd924226e6ff28fc 100644 --- a/src/graph/spectral/graph_nonbacktracking.hh +++ b/src/graph/spectral/graph_nonbacktracking.hh @@ -54,6 +54,117 @@ void get_nonbacktracking(Graph& g, Index index, } } +template +void nbt_matvec(Graph& g, EIndex eindex, V& x, V& ret) +{ + auto get_idx = + [&](auto& e, bool reverse = false) + { + int64_t idx = get(eindex, e); + if constexpr (!graph_tool::is_directed_::apply::type::value) + { + size_t u = source(e, g); + size_t v = target(e, g); + if (reverse) + std::swap(u, v); + if constexpr (!transpose) + idx = (idx << 1) + (u > v); + else + idx = (idx << 1) + (v > u); + } + return idx; + }; + + + parallel_edge_loop + (g, + [&](const auto& e) + { + size_t u = source(e, g); + size_t v = target(e, g); + + int64_t idx = get_idx(e); + + for (const auto& e2 : out_edges_range(v, g)) + { + auto w = target(e2, g); + if (w == u || w == v) + continue; + int64_t idx2 = get_idx(e2); + ret[idx] += x[idx2]; + } + + idx = get_idx(e, true); + + for (const auto& e2 : out_edges_range(u, g)) + { + auto w = target(e2, g); + if (w == u || w == v) + continue; + int64_t idx2 = get_idx(e2); + ret[idx] += x[idx2]; + } + + }); +} + +template +void nbt_matmat(Graph& g, EIndex eindex, V& x, V& ret) +{ + auto get_idx = + [&](auto& e, bool reverse = false) + { + int64_t idx = get(eindex, e); + if constexpr (!graph_tool::is_directed_::apply::type::value) + { + size_t u = source(e, g); + size_t v = target(e, g); + if (reverse) + std::swap(u, v); + if constexpr (!transpose) + idx = (idx << 1) + (u > v); + else + idx = (idx << 1) + (v > u); + } + return idx; + }; + + + size_t k = x.shape()[1]; + parallel_edge_loop + (g, + [&](const auto& e) + { + size_t u = source(e, g); + size_t v = target(e, g); + + int64_t idx = get_idx(e); + + for (const auto& e2 : out_edges_range(v, g)) + { + auto w = target(e2, g); + if (w == u || w == v) + continue; + int64_t idx2 = get_idx(e2); + for (size_t i = 0; i < k; ++i) + ret[idx][i] += x[idx2][i]; + } + + idx = get_idx(e, true); + + for (const auto& e2 : out_edges_range(u, g)) + { + auto w = target(e2, g); + if (w == u || w == v) + continue; + int64_t idx2 = get_idx(e2); + for (size_t i = 0; i < k; ++i) + ret[idx][i] += x[idx2][i]; + } + + }); +} + template void get_compact_nonbacktracking(Graph& g, std::vector& i, @@ -90,6 +201,76 @@ void get_compact_nonbacktracking(Graph& g, } } +template +void cnbt_matvec(Graph& g, VIndex vindex, V& x, V& ret) +{ + size_t N = HardNumVertices()(g); + parallel_vertex_loop + (g, + [&](const auto& v) + { + size_t i = get(vindex, v); + auto& y = ret[i]; + for (const auto& e : out_edges_range(v, g)) + { + auto u = target(e, g); + size_t j = get(vindex, u); + y += x[j]; + } + + if constexpr (!transpose) + { + ret[i] -= x[i + N]; + ret[i + N] = (out_degree(v, g) - 1) * x[i]; + } + else + { + ret[i + N] -= x[i]; + ret[i] = (out_degree(v, g) - 1) * x[i + N]; + } + }); +} + +template +void cnbt_matmat(Graph& g, VIndex vindex, V& x, V& ret) +{ + size_t k = x.shape()[1]; + size_t N = HardNumVertices()(g); + parallel_vertex_loop + (g, + [&](const auto& v) + { + size_t i = get(vindex, v); + auto y = ret[i]; + for (const auto& e : out_edges_range(v, g)) + { + auto u = target(e, g); + size_t j = get(vindex, u); + for (size_t l = 0; l < k; ++l) + y[l] += x[j][l]; + } + + auto d = (out_degree(v, g) - 1); + if constexpr (!transpose) + { + for (size_t l = 0; l < k; ++l) + { + ret[i][l] -= x[i + N][l]; + ret[i + N][l] = d * x[i][l]; + } + } + else + { + for (size_t l = 0; l < k; ++l) + { + ret[i + N][l] -= x[i][l]; + ret[i][l] = d * x[i + N][l]; + } + } + }); +} + + } // namespace graph_tool #endif // GRAPH_NONBACKTRACKING_MATRIX_HH diff --git a/src/graph/spectral/graph_norm_laplacian.cc b/src/graph/spectral/graph_norm_laplacian.cc index aaaf241e6b5d1d6970b08f64dbc4d5ae151d9986..5cc478c9b830a52f1866f38b4f09551c8ef97e12 100644 --- a/src/graph/spectral/graph_norm_laplacian.cc +++ b/src/graph/spectral/graph_norm_laplacian.cc @@ -70,3 +70,67 @@ void norm_laplacian(GraphInterface& g, boost::any index, boost::any weight, }, vertex_scalar_properties(), weight_props_t())(index, weight); } + +void norm_laplacian_matvec(GraphInterface& g, boost::any index, + boost::any weight, boost::any deg, python::object ov, + python::object oret) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(deg).get_unchecked(); + + run_action<>() + (g, + [&](auto&& graph, auto&& vi, auto&& w) + { + return nlap_matvec(graph, vi, w, d, v, ret); + }, + vertex_scalar_properties(), weight_props_t())(index, weight); +} + +void norm_laplacian_matmat(GraphInterface& g, boost::any index, + boost::any weight, boost::any deg, python::object ov, + python::object oret) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(deg).get_unchecked(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& vi, auto&& w) + { + return nlap_matmat(graph, vi, w, d, v, ret); + }, + vertex_scalar_properties(), weight_props_t())(index, weight); +} diff --git a/src/graph/spectral/graph_transition.cc b/src/graph/spectral/graph_transition.cc index 87e6f1e9109274ad701f24e63388833b5dd7c187..a0edae1b0f2f1a44a8b0f9c0ef124514385e13c9 100644 --- a/src/graph/spectral/graph_transition.cc +++ b/src/graph/spectral/graph_transition.cc @@ -61,3 +61,73 @@ void transition(GraphInterface& g, boost::any index, boost::any weight, }, vertex_scalar_properties(), weight_props_t())(index, weight); } + +void transition_matvec(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(deg).get_unchecked(); + + run_action<>() + (g, + [&](auto&& graph, auto&& vi, auto&& w) + { + if (!transpose) + return trans_matvec(graph, vi, w, d, v, ret); + else + return trans_matvec(graph, vi, w, d, v, ret); + }, + vertex_scalar_properties(), weight_props_t())(index, weight); +} + +void transition_matmat(GraphInterface& g, boost::any index, boost::any weight, + boost::any deg, python::object ov, python::object oret, + bool transpose) +{ + if (!belongs()(index)) + throw ValueException("index vertex property must have a scalar value type"); + + typedef UnityPropertyMap weight_map_t; + typedef mpl::push_back::type + weight_props_t; + + if (!weight.empty() && !belongs()(weight)) + throw ValueException("weight edge property must have a scalar value type"); + + if(weight.empty()) + weight = weight_map_t(); + + typedef typename vprop_map_t::type deg_t; + deg_t::unchecked_t d = any_cast(deg).get_unchecked(); + + multi_array_ref v = get_array(ov); + multi_array_ref ret = get_array(oret); + + run_action<>() + (g, + [&](auto&& graph, auto&& vi, auto&& w) + { + if (!transpose) + return trans_matmat(graph, vi, w, d, v, ret); + else + return trans_matmat(graph, vi, w, d, v, ret); + }, + vertex_scalar_properties(), weight_props_t())(index, weight); +} diff --git a/src/graph/spectral/graph_transition.hh b/src/graph/spectral/graph_transition.hh index 5bb23760252c1340cea95f584c48f7aa7fa6d036..7cba3354c10cb9ec3f6d539313e133db0fb57d7b 100644 --- a/src/graph/spectral/graph_transition.hh +++ b/src/graph/spectral/graph_transition.hh @@ -69,6 +69,93 @@ struct get_transition } }; +template +void trans_matvec(Graph& g, Vindex vindex, Weight w, Deg id, V& x, V& ret) +{ + parallel_vertex_loop + (g, + [&](auto v) + { + std::remove_reference_t y = 0; + if constexpr (!std::is_same_v>) + { + for (auto e : in_or_out_edges_range(v, g)) + { + auto u = target(e, g); + auto w_e = get(w, e); + auto ui = get(vindex, u); + if constexpr (!transpose) + y += w_e * x[ui] * id[u]; + else + y += w_e * x[ui]; + } + } + else + { + for (auto u : in_or_out_neighbors_range(v, g)) + { + auto ui = get(vindex, u); + if constexpr (!transpose) + y += x[ui] * id[u]; + else + y += x[ui]; + } + } + if constexpr (transpose) + y *= id[v]; + ret[get(vindex, v)] = y; + }); +} + + +template +void trans_matmat(Graph& g, Vindex vindex, Weight w, Deg id, Mat& x, Mat& ret) +{ + size_t M = x.shape()[1]; + parallel_vertex_loop + (g, + [&](auto v) + { + auto vi = get(vindex, v); + auto y = ret[vi]; + if constexpr (!std::is_same_v>) + { + for (auto e : in_or_out_edges_range(v, g)) + { + auto u = target(e, g); + auto w_e = get(w, e); + auto ui = get(vindex, u); + for (size_t i = 0; i < M; ++i) + { + if constexpr (!transpose) + y[i] += w_e * x[ui][i] * id[u]; + else + y[i] += w_e * x[ui][i]; + } + } + } + else + { + for (auto u : in_or_out_neighbors_range(v, g)) + { + auto ui = get(vindex, u); + for (size_t i = 0; i < M; ++i) + { + if constexpr (!transpose) + y[i] += x[ui][i] * id[u]; + else + y[i] += x[ui][i]; + } + } + } + if constexpr (transpose) + { + for (size_t i = 0; i < M; ++i) + y[i] *= id[v]; + } + }); +} + } // namespace graph_tool #endif // GRAPH_TRANSITION_HH diff --git a/src/graph_tool/spectral/__init__.py b/src/graph_tool/spectral/__init__.py index 5bdda5685ef61c4a73f21d0a18905acd53e15563..010d01f5aedf36f45d61001b567482c711539345 100644 --- a/src/graph_tool/spectral/__init__.py +++ b/src/graph_tool/spectral/__init__.py @@ -35,6 +35,13 @@ Summary modularity_matrix hashimoto + AdjacencyOperator + LaplacianOperator + IncidenceOperator + TransitionOperator + HashimotoOperator + CompactHashimotoOperator + Contents ++++++++ """ @@ -49,11 +56,12 @@ import scipy.sparse.linalg from .. dl_import import dl_import dl_import("from . import libgraph_tool_spectral") -__all__ = ["adjacency", "laplacian", "incidence", "transition", - "modularity_matrix", "hashimoto"] - +__all__ = ["adjacency", "AdjacencyOperator", "laplacian", "LaplacianOperator", + "incidence", "IncidenceOperator", "transition", "TransitionOperator", + "modularity_matrix", "hashimoto", "HashimotoOperator", + "CompactHashimotoOperator"] -def adjacency(g, weight=None, index=None): +def adjacency(g, weight=None, vindex=None, operator=False): r"""Return the adjacency matrix of the graph. Parameters @@ -62,13 +70,16 @@ def adjacency(g, weight=None, index=None): Graph to be used. weight : :class:`~graph_tool.EdgePropertyMap` (optional, default: True) Edge property map with the edge weights. - index : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) + vindex : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) Vertex property map specifying the row/column indexes. If not provided, the internal vertex index is used. + operator : ``bool`` (optional, default: ``False``) + If ``True``, a :class:`scipy.sparse.linalg.LinearOperator` subclass is + returned, instead of a sparse matrix. Returns ------- - A : :class:`~scipy.sparse.csr_matrix` + A : :class:`~scipy.sparse.csr_matrix` or :class:`AdjacencyOperator` The (sparse) adjacency matrix. Notes @@ -110,6 +121,17 @@ def adjacency(g, weight=None, index=None): i`. Although this is a typical definition in network and graph theory literature, many also use the transpose of this matrix. + .. note:: + + For many linear algebra computations it is more efficient to pass + ``operator=True``. This makes this function return a + :class:`scipy.sparse.linalg.LinearOperator` subclass, which implements + matrix-vector and matrix-matrix multiplication, and is sufficient for + the sparse linear algebra operations available in the scipy module + :mod:`scipy.sparse.linalg`. This avoids copying the whole graph as a + sparse matrix, and performs the multiplication operations in parallel + (if enabled during compilation). + Examples -------- .. testsetup:: @@ -118,8 +140,11 @@ def adjacency(g, weight=None, index=None): from pylab import * >>> g = gt.collection.data["polblogs"] - >>> A = gt.adjacency(g) - >>> ew, ev = scipy.linalg.eig(A.todense()) + >>> A = gt.adjacency(g, operator=True) + >>> N = g.num_vertices() + >>> ew1 = scipy.sparse.linalg.eigs(A, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(A, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) >>> figure(figsize=(8, 2)) <...> @@ -135,7 +160,7 @@ def adjacency(g, weight=None, index=None): .. figure:: adjacency-spectrum.* :align: center - Adjacency matrix spectrum for the political blog network. + Adjacency matrix spectrum for the political blogs network. References ---------- @@ -143,12 +168,15 @@ def adjacency(g, weight=None, index=None): """ - if index is None: + if operator: + return AdjacencyOperator(g, weight=weight, vindex=vindex) + + if vindex is None: if g.get_vertex_filter()[0] is not None: - index = g.new_vertex_property("int64_t") - index.fa = numpy.arange(g.num_vertices()) + vindex = g.new_vertex_property("int64_t") + vindex.fa = numpy.arange(g.num_vertices()) else: - index = g.vertex_index + vindex = g.vertex_index E = g.num_edges() if g.is_directed() else 2 * g.num_edges() @@ -156,7 +184,7 @@ def adjacency(g, weight=None, index=None): i = numpy.zeros(E, dtype="int32") j = numpy.zeros(E, dtype="int32") - libgraph_tool_spectral.adjacency(g._Graph__graph, _prop("v", g, index), + libgraph_tool_spectral.adjacency(g._Graph__graph, _prop("v", g, vindex), _prop("e", g, weight), data, i, j) if E > 0: @@ -168,9 +196,58 @@ def adjacency(g, weight=None, index=None): m = m.tocsr() return m +class AdjacencyOperator(scipy.sparse.linalg.LinearOperator): + def __init__(self, g, weight=None, vindex=None): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the adjacency + matrix of a graph. See :func:`adjacency` for details.""" + self.g = g + self.weight = weight + + if vindex is None: + if g.get_vertex_filter()[0] is not None: + self.vindex = g.new_vertex_property("int64_t") + self.vindex.fa = numpy.arange(g.num_vertices()) + N = g.num_vertices() + else: + self.vindex = g.vertex_index + N = g.num_vertices() + else: + self.vindex = vindex + if vindex is vindex.get_graph().vertex_index: + N = g.num_vertices() + else: + N = vindex.fa.max() + 1 + + self.shape = (N, N) + self.dtype = numpy.dtype("float") + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.adjacency_matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + x, y) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.adjacency_matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + x, y) + return y + + def _adjoint(self): + if self.g.is_directed(): + return AdjacencyOperator(GraphView(self.g, reversed=True), + self.weight, self.vindex) + else: + return self @_limit_args({"deg": ["total", "in", "out"]}) -def laplacian(g, deg="total", normalized=False, weight=None, index=None): +def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False): r"""Return the Laplacian matrix of the graph. Parameters @@ -179,17 +256,20 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): Graph to be used. deg : str (optional, default: "total") Degree to be used, in case of a directed graph. - normalized : bool (optional, default: False) + norm : bool (optional, default: False) Whether to compute the normalized Laplacian. weight : :class:`~graph_tool.EdgePropertyMap` (optional, default: True) Edge property map with the edge weights. - index : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) + vindex : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) Vertex property map specifying the row/column indexes. If not provided, the internal vertex index is used. + operator : ``bool`` (optional, default: ``False``) + If ``True``, a :class:`scipy.sparse.linalg.LinearOperator` subclass is + returned, instead of a sparse matrix. Returns ------- - l : :class:`~scipy.sparse.csr_matrix` + L : :class:`~scipy.sparse.csr_matrix` or :class:`LaplacianOperator` The (sparse) Laplacian matrix. Notes @@ -221,7 +301,7 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): 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_{ji}w_{ji}` - if ``deg=="out"`` or :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij}` ``deg=="in"``. + if ``deg=="out"`` or :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij}` if ``deg=="in"``. .. note:: @@ -230,6 +310,17 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): i`. Although this is a typical definition in network and graph theory literature, many also use the transpose of this matrix. + .. note:: + + For many linear algebra computations it is more efficient to pass + ``operator=True``. This makes this function return a + :class:`scipy.sparse.linalg.LinearOperator` subclass, which implements + matrix-vector and matrix-matrix multiplication, and is sufficient for + the sparse linear algebra operations available in the scipy module + :mod:`scipy.sparse.linalg`. This avoids copying the whole graph as a + sparse matrix, and performs the multiplication operations in parallel + (if enabled during compilation). + Examples -------- @@ -239,8 +330,11 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): from pylab import * >>> g = gt.collection.data["polblogs"] - >>> L = gt.laplacian(g) - >>> ew, ev = scipy.linalg.eig(L.todense()) + >>> L = gt.laplacian(g, operator=True) + >>> N = g.num_vertices() + >>> ew1 = scipy.sparse.linalg.eigs(L, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(L, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) >>> figure(figsize=(8, 2)) <...> @@ -256,10 +350,12 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): .. figure:: laplacian-spectrum.* :align: center - Laplacian matrix spectrum for the political blog network. + Laplacian matrix spectrum for the political blogs network. - >>> L = gt.laplacian(g, normalized=True) - >>> ew, ev = scipy.linalg.eig(L.todense()) + >>> L = gt.laplacian(g, norm=True, operator=True) + >>> ew1 = scipy.sparse.linalg.eigs(L, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(L, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) >>> figure(figsize=(8, 2)) <...> @@ -275,19 +371,23 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): .. figure:: norm-laplacian-spectrum.* :align: center - Normalized Laplacian matrix spectrum for the political blog network. + Normalized Laplacian matrix spectrum for the political blogs network. References ---------- .. [wikipedia-laplacian] http://en.wikipedia.org/wiki/Laplacian_matrix """ - if index is None: + if operator: + return LaplacianOperator(g, deg=deg, norm=norm, weight=weight, + vindex=vindex) + + if vindex is None: if g.get_vertex_filter()[0] is not None: - index = g.new_vertex_property("int64_t") - index.fa = numpy.arange(g.num_vertices()) + vindex = g.new_vertex_property("int64_t") + vindex.fa = numpy.arange(g.num_vertices()) else: - index = g.vertex_index + vindex = g.vertex_index V = g.num_vertices() nself = int(label_self_loops(g, mark_only=True).a.sum()) @@ -300,11 +400,11 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): 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), + if norm: + libgraph_tool_spectral.norm_laplacian(g._Graph__graph, _prop("v", g, vindex), _prop("e", g, weight), deg, data, i, j) else: - libgraph_tool_spectral.laplacian(g._Graph__graph, _prop("v", g, index), + libgraph_tool_spectral.laplacian(g._Graph__graph, _prop("v", g, vindex), _prop("e", g, weight), deg, data, i, j) if E > 0: V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1)) @@ -315,8 +415,85 @@ def laplacian(g, deg="total", normalized=False, weight=None, index=None): m = m.tocsr() return m +class LaplacianOperator(scipy.sparse.linalg.LinearOperator): + + @_limit_args({"deg": ["total", "in", "out"]}) + def __init__(self, g, weight=None, deg="out", norm=False, vindex=None): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the laplacian + matrix of a graph. See :func:`laplacian` for details.""" + + self.g = g + self.weight = weight -def incidence(g, vindex=None, eindex=None): + if vindex is None: + if g.get_vertex_filter()[0] is not None: + self.vindex = g.new_vertex_property("int64_t") + self.vindex.fa = numpy.arange(g.num_vertices()) + N = g.num_vertices() + else: + self.vindex = g.vertex_index + N = g.num_vertices() + else: + self.vindex = vindex + if vindex is vindex.get_graph().vertex_index: + N = g.num_vertices() + else: + N = vindex.fa.max() + 1 + + self.shape = (N, N) + self.deg = deg + self.d = self.g.degree_property_map(deg, weight) + if norm: + idx = self.d.a > 0 + d = g.new_vp("double") + d.a[idx] = 1./numpy.sqrt(self.d.a[idx]) + self.d = d + else: + self.d = self.d.copy("double") + self.norm = norm + self.dtype = numpy.dtype("float") + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + if self.norm: + matvec = libgraph_tool_spectral.norm_laplacian_matvec + else: + matvec = libgraph_tool_spectral.laplacian_matvec + matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + _prop("v", self.g, self.d), x, y) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + if self.norm: + matmat = libgraph_tool_spectral.norm_laplacian_matmat + else: + matmat = libgraph_tool_spectral.laplacian_matmat + matmat(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + _prop("v", self.g, self.d), x, y) + return y + + def _adjoint(self): + if self.g.is_directed(): + deg = self.deg + if deg == "in": + deg = "out" + elif deg == "out": + deg = "in" + return LaplacianOperator(GraphView(self.g, reversed=True), + self.weight, deg=deg, norm=self.norm, + vindex=self.vindex) + else: + return self + + +def incidence(g, vindex=None, eindex=None, operator=False): r"""Return the incidence matrix of the graph. Parameters @@ -329,10 +506,13 @@ def incidence(g, vindex=None, eindex=None): eindex : :class:`~graph_tool.EdgePropertyMap` (optional, default: None) Edge property map specifying the column indexes. If not provided, the internal edge index is used. + operator : ``bool`` (optional, default: ``False``) + If ``True``, a :class:`scipy.sparse.linalg.LinearOperator` subclass is + returned, instead of a sparse matrix. Returns ------- - a : :class:`~scipy.sparse.csr_matrix` + a : :class:`~scipy.sparse.csr_matrix` or :class:`IncidenceOperator` The (sparse) incidence matrix. Notes @@ -358,28 +538,57 @@ def incidence(g, vindex=None, eindex=None): 0 & \text{otherwise} \end{cases} + .. note:: + + For many linear algebra computations it is more efficient to pass + ``operator=True``. This makes this function return a + :class:`scipy.sparse.linalg.LinearOperator` subclass, which implements + matrix-vector and matrix-matrix multiplication, and is sufficient for + the sparse linear algebra operations available in the scipy module + :mod:`scipy.sparse.linalg`. This avoids copying the whole graph as a + sparse matrix, and performs the multiplication operations in parallel + (if enabled during compilation). + Examples -------- + .. testsetup:: - gt.seed_rng(42) + import scipy.linalg + from pylab import * + + >>> g = gt.collection.data["polblogs"] + >>> B = gt.incidence(g, operator=True) + >>> N = g.num_vertices() + >>> s1 = scipy.sparse.linalg.svds(B, k=N//2, which="LM", return_singular_vectors=False) + >>> s2 = scipy.sparse.linalg.svds(B, k=N-N//2, which="SM", return_singular_vectors=False) + >>> s = np.concatenate((s1, s2)) + >>> s.sort() + + >>> figure(figsize=(8, 2)) + <...> + >>> plot(s, "s") + [...] + >>> xlabel(r"$i$") + Text(...) + >>> ylabel(r"$\lambda_i$") + Text(...) + >>> tight_layout() + >>> savefig("polblogs-indidence-svd.svg") + + .. figure:: polblogs-indidence-svd.* + :align: center - >>> g = gt.random_graph(100, lambda: (2,2)) - >>> m = gt.incidence(g) - >>> print(m.todense()) - [[-1. -1. 0. ... 0. 0. 0.] - [ 0. 0. 0. ... 0. 0. 0.] - [ 1. 0. 0. ... 0. 0. 0.] - ... - [ 0. 0. -1. ... 0. 0. 0.] - [ 0. 0. 0. ... 0. 0. 0.] - [ 0. 0. 0. ... 0. 0. 0.]] + Incidence singular values for the political blogs network. References ---------- .. [wikipedia-incidence] http://en.wikipedia.org/wiki/Incidence_matrix """ + if operator: + return IncidenceOperator(g, vindex=vindex, eindex=eindex) + if vindex is None: if g.get_edge_filter()[0] is not None: vindex = g.new_vertex_property("int64_t") @@ -409,7 +618,76 @@ def incidence(g, vindex=None, eindex=None): m = m.tocsr() return m -def transition(g, weight=None, index=None): + +class IncidenceOperator(scipy.sparse.linalg.LinearOperator): + + def __init__(self, g, vindex=None, eindex=None, transpose=False): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the incidence + matrix of a graph. See :func:`incidence` for details. + """ + + self.g = g + self.transpose = transpose + self.dtype = numpy.dtype("float") + + if vindex is None: + if g.get_vertex_filter()[0] is not None: + self.vindex = g.new_vertex_property("int64_t") + self.vindex.fa = numpy.arange(g.num_vertices()) + N = g.num_vertices() + else: + self.vindex = g.vertex_index + N = g.num_vertices() + else: + self.vindex = vindex + if vindex is vindex.get_graph().vertex_index: + N = g.num_vertices() + else: + N = vindex.fa.max() + 1 + + if eindex is None: + if g.get_edge_filter()[0] is not None: + self.eindex = g.new_edge_property("int64_t") + self.eindex.fa = numpy.arange(g.num_edges()) + E = g.num_edges() + else: + self.eindex = g.edge_index + E = g.edge_index_range + else: + self.eindex = eindex + if eindex is g.edge_index: + E = g.edge_index_range + else: + E = self.eindex.fa.max() + 1 + + if not transpose: + self.shape = (N, E) + else: + self.shape = (E, N) + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.incidence_matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.eindex), + x, y, self.transpose) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.incidence_matmat(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.eindex), + x, y, self.transpose) + return y + + def _adjoint(self): + return IncidenceOperator(self.g, vindex=self.vindex, eindex=self.eindex, + transpose=not self.transpose) + +def transition(g, weight=None, vindex=None, operator=False): r"""Return the transition matrix of the graph. Parameters @@ -418,13 +696,16 @@ def transition(g, weight=None, index=None): Graph to be used. weight : :class:`~graph_tool.EdgePropertyMap` (optional, default: True) Edge property map with the edge weights. - index : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) - Vertex property map specifying the row/column indexes. If not provided, the - internal vertex index is used. + vindex : :class:`~graph_tool.VertexPropertyMap` (optional, default: None) + Vertex property map specifying the row/column indexes. If not provided, + the internal vertex index is used. + operator : ``bool`` (optional, default: ``False``) + If ``True``, a :class:`scipy.sparse.linalg.LinearOperator` subclass is + returned, instead of a sparse matrix. Returns ------- - T : :class:`~scipy.sparse.csr_matrix` + T : :class:`~scipy.sparse.csr_matrix` or :class:`TransitionOperator` The (sparse) transition matrix. Notes @@ -448,6 +729,16 @@ def transition(g, weight=None, index=None): i`. Although this is a typical definition in network and graph theory literature, many also use the transpose of this matrix. + .. note:: + + For many linear algebra computations it is more efficient to pass + ``operator=True``. This makes this function return a + :class:`scipy.sparse.linalg.LinearOperator` subclass, which implements + matrix-vector and matrix-matrix multiplication, and is sufficient for + the sparse linear algebra operations available in the scipy module + :mod:`scipy.sparse.linalg`. This avoids copying the whole graph as a + sparse matrix, and performs the multiplication operations in parallel + (if enabled during compilation). Examples -------- @@ -457,8 +748,11 @@ def transition(g, weight=None, index=None): from pylab import * >>> g = gt.collection.data["polblogs"] - >>> T = gt.transition(g) - >>> ew, ev = scipy.linalg.eig(T.todense()) + >>> T = gt.transition(g, operator=True) + >>> N = g.num_vertices() + >>> ew1 = scipy.sparse.linalg.eigs(T, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(T, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) >>> figure(figsize=(8, 2)) <...> @@ -474,26 +768,30 @@ def transition(g, weight=None, index=None): .. figure:: transition-spectrum.* :align: center - Transition matrix spectrum for the political blog network. + Transition matrix spectrum for the political blogs network. References ---------- .. [wikipedia-transition] https://en.wikipedia.org/wiki/Stochastic_matrix + """ - if index is None: + if operator: + return TransitionOperator(g, weight=weight, vindex=vindex) + + if vindex is None: if g.get_vertex_filter()[0] is not None: - index = g.new_vertex_property("int64_t") - index.fa = numpy.arange(g.num_vertices()) + vindex = g.new_vertex_property("int64_t") + vindex.fa = numpy.arange(g.num_vertices()) else: - index = g.vertex_index + vindex = 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), + libgraph_tool_spectral.transition(g._Graph__graph, _prop("v", g, vindex), _prop("e", g, weight), data, i, j) if E > 0: @@ -504,9 +802,74 @@ def transition(g, weight=None, index=None): m = m.tocsr() return m +class TransitionOperator(scipy.sparse.linalg.LinearOperator): + + def __init__(self, g, weight=None, transpose=False, inv_d=None, vindex=None): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the transition + matrix of a graph. See :func:`transition` for details. + """ + + self.g = g + self.weight = weight + + if vindex is None: + if g.get_vertex_filter()[0] is not None: + self.vindex = g.new_vertex_property("int64_t") + self.vindex.fa = numpy.arange(g.num_vertices()) + N = g.num_vertices() + else: + self.vindex = g.vertex_index + N = g.num_vertices() + else: + self.vindex = vindex + if vindex is vindex.get_graph().vertex_index: + N = g.num_vertices() + else: + N = vindex.fa.max() + 1 + + self.shape = (N, N) + if inv_d is None: + d = self.g.degree_property_map("out", weight) + nd = g.new_vp("double") + idx = d.a > 0 + nd.a[idx] = 1/d.a[idx] + self.d = nd + else: + self.d = inv_d.copy("double") + self.dtype = numpy.dtype("float") + self.transpose = transpose + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.transition_matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + _prop("v", self.g, self.d), x, + y, self.transpose) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.transition_matmat(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + _prop("e", self.g, self.weight), + _prop("v", self.g, self.d), x, + y, self.transpose) + return y + + def _adjoint(self): + if self.g.is_directed(): + g = GraphView(self.g, reversed=True) + else: + g = self.g + return TransitionOperator(g, self.weight, inv_d=self.d, + transpose=not self.transpose, + vindex=self.vindex) -def modularity_matrix(g, weight=None, index=None): +def modularity_matrix(g, weight=None, vindex=None): r"""Return the modularity matrix of the graph. Parameters @@ -549,8 +912,10 @@ def modularity_matrix(g, weight=None, index=None): >>> 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) + >>> N = g.num_vertices() + >>> ew1 = scipy.sparse.linalg.eigs(B, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(B, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) >>> figure(figsize=(8, 2)) <...> @@ -560,13 +925,14 @@ def modularity_matrix(g, weight=None, index=None): Text(...) >>> ylabel(r"$\operatorname{Im}(\lambda)$") Text(...) + >>> autoscale() >>> tight_layout() >>> savefig("modularity-spectrum.svg") .. figure:: modularity-spectrum.* :align: center - Modularity matrix spectrum for the political blog network. + Modularity matrix spectrum for the political blogs network. References ---------- @@ -575,28 +941,29 @@ def modularity_matrix(g, weight=None, index=None): :doi:`10.1103/PhysRevE.69.026113` """ - A = adjacency(g, weight=weight, index=index) + A = adjacency(g, weight=weight, vindex=vindex, operator=True) + A_T = A.adjoint() 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] + N = g.num_vertices() 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 + nx = A.matvec(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 + nx = A_T.matvec(x) - k_in * numpy.dot(k_out, x) / E2 return nx B = scipy.sparse.linalg.LinearOperator((g.num_vertices(), g.num_vertices()), @@ -605,7 +972,7 @@ def modularity_matrix(g, weight=None, index=None): return B -def hashimoto(g, index=None, compact=False): +def hashimoto(g, index=None, compact=False, operator=False): r"""Return the Hashimoto (or non-backtracking) matrix of a graph. Parameters @@ -618,10 +985,13 @@ def hashimoto(g, index=None, compact=False): compact : ``boolean`` (optional, default: ``False``) If ``True``, a compact :math:`2|V|\times 2|V|` version of the matrix is returned. + operator : ``bool`` (optional, default: ``False``) + If ``True``, a :class:`scipy.sparse.linalg.LinearOperator` subclass is + returned, instead of a sparse matrix. Returns ------- - H : :class:`~scipy.sparse.csr_matrix` + H : :class:`~scipy.sparse.csr_matrix` or :class:`HashimotoOperator` or :class:`CompactHashimotoOperator` The (sparse) Hashimoto matrix. Notes @@ -655,6 +1025,17 @@ def hashimoto(g, index=None, compact=False): where :math:`\boldsymbol A` is the adjacency matrix, and :math:`\boldsymbol D` is the diagonal matrix with the node degrees [krzakala_spectral]_. + .. note:: + + For many linear algebra computations it is more efficient to pass + ``operator=True``. This makes this function return a + :class:`scipy.sparse.linalg.LinearOperator` subclass, which implements + matrix-vector and matrix-matrix multiplication, and is sufficient for + the sparse linear algebra operations available in the scipy module + :mod:`scipy.sparse.linalg`. This avoids copying the whole graph as a + sparse matrix, and performs the multiplication operations in parallel + (if enabled during compilation). + Examples -------- .. testsetup:: @@ -663,8 +1044,12 @@ def hashimoto(g, index=None, compact=False): from pylab import * >>> g = gt.collection.data["football"] - >>> H = gt.hashimoto(g) - >>> ew, ev = scipy.linalg.eig(H.todense()) + >>> H = gt.hashimoto(g, operator=True) + >>> N = 2 * g.num_edges() + >>> ew1 = scipy.sparse.linalg.eigs(H, k=N//2, which="LR", return_eigenvectors=False) + >>> ew2 = scipy.sparse.linalg.eigs(H, k=N-N//2, which="SR", return_eigenvectors=False) + >>> ew = np.concatenate((ew1, ew2)) + >>> figure(figsize=(8, 4)) <...> >>> scatter(real(ew), imag(ew), c=sqrt(abs(ew)), linewidths=0, alpha=0.6) @@ -694,6 +1079,9 @@ def hashimoto(g, index=None, compact=False): """ if compact: + if operator: + return CompactHashimotoOperator(g) + i = Vector_int64_t() j = Vector_int64_t() x = Vector_double() @@ -704,6 +1092,9 @@ def hashimoto(g, index=None, compact=False): N = g.num_vertices(ignore_filter=True) m = scipy.sparse.coo_matrix((x, (i.a,j.a)), shape=(2 * N, 2 * N)) else: + if operator: + return HashimotoOperator(g, eindex=index) + if index is None: if g.get_edge_filter()[0] is not None: index = g.new_edge_property("int64_t") @@ -726,3 +1117,107 @@ def hashimoto(g, index=None, compact=False): m = scipy.sparse.coo_matrix((data, (i.a,j.a)), shape=(E, E)) m = m.tocsr() return m + + +class HashimotoOperator(scipy.sparse.linalg.LinearOperator): + + def __init__(self, g, eindex=None, transpose=False): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the hashimoto + matrix of a graph. See :func:`hashimoto` for details. + """ + + self.g = g + if eindex is None: + if g.get_edge_filter()[0] is not None: + self.eindex = g.new_edge_property("int64_t") + self.eindex.fa = numpy.arange(g.num_edges()) + E = g.num_edges() + else: + self.eindex = g.edge_index + E = g.edge_index_range + else: + self.eindex = eindex + if eindex is g.edge_index: + E = g.edge_index_range + else: + E = self.eindex.fa.max() + 1 + if g.is_directed(): + self.shape = (E, E) + else: + self.shape = (2 * E, 2 * E) + + self.dtype = numpy.dtype("float") + self.transpose = transpose + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.nonbacktracking_matvec(self.g._Graph__graph, + _prop("e", self.g, self.eindex), + x, y, self.transpose) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.nonbacktracking_matmat(self.g._Graph__graph, + _prop("e", self.g, self.eindex), + x, y, self.transpose) + return y + + def _adjoint(self): + if self.g.is_directed(): + g = GraphView(self.g, reversed=True) + else: + g = self.g + return HashimotoOperator(g, self.eindex, transpose=not self.transpose) + +class CompactHashimotoOperator(scipy.sparse.linalg.LinearOperator): + + def __init__(self, g, vindex=None, transpose=False): + r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the compact + hashimoto matrix of a graph. See :func:`hashimoto` for details. + """ + + self.g = g + if vindex is None: + if g.get_vertex_filter()[0] is not None: + self.vindex = g.new_vertex_property("int64_t") + self.vindex.fa = numpy.arange(g.num_vertices()) + N = g.num_vertices() + else: + self.vindex = g.vertex_index + N = g.num_vertices() + else: + self.vindex = vindex + if vindex is vindex.get_graph().vertex_index: + N = g.num_vertices() + else: + N = vindex.fa.max() + 1 + + self.shape = (2 * N, 2 * N) + self.dtype = numpy.dtype("float") + self.transpose = transpose + + def _matvec(self, x): + y = numpy.zeros(self.shape[0]) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.compact_nonbacktracking_matvec(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + x, y, self.transpose) + return y + + def _matmat(self, x): + y = numpy.zeros((self.shape[0], x.shape[1])) + x = numpy.asarray(x, dtype="float") + libgraph_tool_spectral.compact_nonbacktracking_matmat(self.g._Graph__graph, + _prop("v", self.g, self.vindex), + x, y, self.transpose) + return y + + def _adjoint(self): + if self.g.is_directed(): + g = GraphView(self.g, reversed=True) + else: + g = self.g + return CompactHashimotoOperator(g, self.vindex, transpose=not self.transpose)