From cfd1a4e1bf2c40420333c5aca1efea689a0b319d Mon Sep 17 00:00:00 2001 From: Thomas Date: Wed, 20 May 2026 16:20:36 +0200 Subject: [PATCH] feat: add h_sparse_triplets method and corresponding tests for sparse Hessian --- include/hypergraph/hypergraph.h | 113 +++++++++++++++++- pyproject.toml | 1 + src/Module.cpp | 1 + tests/TestHyperGraph.py | 200 ++++++++++++++++++++++++++++++++ 4 files changed, 314 insertions(+), 1 deletion(-) diff --git a/include/hypergraph/hypergraph.h b/include/hypergraph/hypergraph.h index e3422a1..1a7f8d9 100644 --- a/include/hypergraph/hypergraph.h +++ b/include/hypergraph/hypergraph.h @@ -20,6 +20,7 @@ #include #include #include +#include namespace hypergraph { @@ -654,6 +655,84 @@ class HyperGraph { } } + index num_variables() const + { + return length(m_variables); + } + + std::tuple, std::vector, std::vector> + h_sparse_triplets(const bool full = false) const + { + const index n = length(m_variables); + + std::vector rows; + std::vector cols; + std::vector values; + + // Build a map from vertex_id → variable_index for efficient lookup + tsl::robin_map vid_to_var; + vid_to_var.reserve(n); + for (index i = 0; i < n; i++) { + vid_to_var[m_variables[i].id()] = i; + } + + // Diagonal entries from m_self_second_order_edges + for (index i = 0; i < n; i++) { + const T val = m_self_second_order_edges[m_variables[i].id()]; + if (val != T{}) { + rows.push_back(i); + cols.push_back(i); + values.push_back(double(val)); + } + } + + // Off-diagonal entries from m_second_order_edges + // The structure stores entries with key = min(vid_a, vid_b) inside + // m_second_order_edges[max(vid_a, vid_b)]. + for (index outer_vid = 0; outer_vid < length(m_second_order_edges); outer_vid++) { + const auto& btree = m_second_order_edges[outer_vid]; + if (btree.empty()) { + continue; + } + + // Check if the outer vertex_id corresponds to a variable + auto outer_it = vid_to_var.find(outer_vid); + if (outer_it == vid_to_var.end()) { + continue; + } + const index outer_var_idx = outer_it->second; + + for (const auto& [inner_vid, val] : btree) { + if (val == T{}) { + continue; + } + + auto inner_it = vid_to_var.find(inner_vid); + if (inner_it == vid_to_var.end()) { + continue; + } + const index inner_var_idx = inner_it->second; + + // inner_vid < outer_vid by construction (minmax in second_order_edge) + // So inner_var_idx may or may not be < outer_var_idx + const auto [row, col] = hypergraph::minmax(inner_var_idx, outer_var_idx); + + // Upper triangle: row <= col + rows.push_back(row); + cols.push_back(col); + values.push_back(double(val)); + + if (full && row != col) { + rows.push_back(col); + cols.push_back(row); + values.push_back(double(val)); + } + } + } + + return {std::move(rows), std::move(cols), std::move(values)}; + } + public: // python template static void register_python(TModule& m) @@ -666,6 +745,8 @@ class HyperGraph { nb::class_(m, name.c_str()) // constructors .def(nb::init<>()) + // properties + .def_prop_ro("num_variables", &Type::num_variables) // methods .def("new_variable", &Type::new_variable, "value"_a) .def("new_variables", &Type::new_variables, "values"_a) @@ -673,7 +754,37 @@ class HyperGraph { .def("g", nb::overload_cast<>(&Type::g, nb::const_)) .def("g", nb::overload_cast>(&Type::g, nb::const_), "out"_a) .def("h", nb::overload_cast(&Type::h, nb::const_), "full"_a = false) - .def("h", nb::overload_cast, const bool>(&Type::h, nb::const_), "out"_a, "full"_a = false); + .def("h", nb::overload_cast, const bool>(&Type::h, nb::const_), "out"_a, "full"_a = false) + .def("h_sparse_triplets", &Type::h_sparse_triplets, "full"_a = false) + .def("h_sparse", [](const Type& self, const std::string& format, const bool full) -> nb::object { + auto [rows, cols, values] = self.h_sparse_triplets(full); + const auto n = self.num_variables(); + + nb::module_ scipy_sparse = nb::module_::import_("scipy.sparse"); + + // Build numpy arrays from the triplet vectors + nb::module_ np = nb::module_::import_("numpy"); + nb::object np_array = np.attr("array"); + + nb::object py_rows = np_array(nb::cast(rows), "dtype"_a = np.attr("int64")); + nb::object py_cols = np_array(nb::cast(cols), "dtype"_a = np.attr("int64")); + nb::object py_vals = np_array(nb::cast(values), "dtype"_a = np.attr("float64")); + + nb::tuple data_ij = nb::make_tuple(py_vals, nb::make_tuple(py_rows, py_cols)); + nb::tuple shape = nb::make_tuple(n, n); + + nb::object coo = scipy_sparse.attr("coo_matrix")(data_ij, "shape"_a = shape); + + if (format == "coo") { + return coo; + } else if (format == "csc") { + return coo.attr("tocsc")(); + } else if (format == "csr") { + return coo.attr("tocsr")(); + } else { + throw std::invalid_argument("h_sparse: format must be 'coo', 'csc', or 'csr'"); + } + }, "format"_a = "csc", "full"_a = false); } }; diff --git a/pyproject.toml b/pyproject.toml index a3eab33..6e733e9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ classifiers = [ ] dependencies = [ "numpy", + "scipy", ] [project.urls] diff --git a/src/Module.cpp b/src/Module.cpp index 914204f..7a3e631 100644 --- a/src/Module.cpp +++ b/src/Module.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include diff --git a/tests/TestHyperGraph.py b/tests/TestHyperGraph.py index 9e86607..c17e3c6 100644 --- a/tests/TestHyperGraph.py +++ b/tests/TestHyperGraph.py @@ -2,6 +2,7 @@ import math import hypergraph as hg import numpy as np +import scipy.sparse from numpy.testing import assert_equal, assert_array_equal, assert_almost_equal, assert_array_almost_equal @@ -1154,5 +1155,204 @@ def test_hypot_pythagorean(self): assert_almost_equal(result.value, 13) + # sparse Hessian + + def test_h_sparse_triplets_basic(self): + """h_sparse_triplets returns correct (rows, cols, values) for upper triangle""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + + result = a * b + graph.compute(result) + + rows, cols, values = graph.h_sparse_triplets() + # a*b has H = [[0, 1], [0, 0]] (upper triangle only) + # Only one non-zero entry: (0, 1) = 1 + assert_equal(len(rows), 1) + assert_equal(len(cols), 1) + assert_equal(len(values), 1) + assert_equal(rows[0], 0) + assert_equal(cols[0], 1) + assert_almost_equal(values[0], 1.0) + + def test_h_sparse_triplets_full(self): + """h_sparse_triplets with full=True returns both triangles""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + + result = a * b + graph.compute(result) + + rows, cols, values = graph.h_sparse_triplets(full=True) + # full=True: (0,1)=1 and (1,0)=1 + assert_equal(len(rows), 2) + # Convert to a set of (row, col, val) for order-independent check + entries = set(zip(rows, cols, values)) + self.assertIn((0, 1, 1.0), entries) + self.assertIn((1, 0, 1.0), entries) + + def test_h_sparse_triplets_diagonal(self): + """h_sparse_triplets includes diagonal entries""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([3, 5]) + + result = a.square() + graph.compute(result) + + rows, cols, values = graph.h_sparse_triplets() + # d²(x²)/dx² = 2, only diagonal entry (0,0) + assert_equal(len(rows), 1) + assert_equal(rows[0], 0) + assert_equal(cols[0], 0) + assert_almost_equal(values[0], 2.0) + + def test_h_sparse_triplets_empty(self): + """h_sparse_triplets returns empty lists when Hessian is zero""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + + result = a + b + graph.compute(result) + + rows, cols, values = graph.h_sparse_triplets() + assert_equal(len(rows), 0) + assert_equal(len(cols), 0) + assert_equal(len(values), 0) + + def test_h_sparse_matches_dense(self): + """h_sparse().toarray() matches h(full=True) for a complex expression""" + graph = hg.HyperGraph() + + ax, ay, az, bx, by, bz = graph.new_variables([1, 2, 3, 4, 5, 6]) + + a = np.array([ax, ay, az]) + b = np.array([bx, by, bz]) + + result = np.linalg.norm(np.cross(a, b)) + + graph.compute(result) + + h_dense = graph.h(full=True) + h_sparse = graph.h_sparse(full=True) + + assert_array_almost_equal(h_sparse.toarray(), h_dense) + + def test_h_sparse_upper_triangle_matches(self): + """h_sparse() default (upper triangle) matches h() upper triangle""" + graph = hg.HyperGraph() + + ax, ay, az, bx, by, bz = graph.new_variables([1, 2, 3, 4, 5, 6]) + + a = np.array([ax, ay, az]) + b = np.array([bx, by, bz]) + + result = np.linalg.norm(np.cross(a, b)) + + graph.compute(result) + + h_dense_upper = graph.h(full=False) + h_sparse_upper = graph.h_sparse(full=False) + + # The sparse version should have the upper triangle entries + # When converted to dense, lower triangle should be zero + sparse_dense = h_sparse_upper.toarray() + assert_array_almost_equal(sparse_dense, h_dense_upper) + + def test_h_sparse_coo_format(self): + """h_sparse(format='coo') returns COO matrix""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + h_sparse = graph.h_sparse(format='coo') + self.assertIsInstance(h_sparse, scipy.sparse.coo_matrix) + + def test_h_sparse_csc_format(self): + """h_sparse(format='csc') returns CSC matrix (default)""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + h_sparse = graph.h_sparse(format='csc') + self.assertIsInstance(h_sparse, scipy.sparse.csc_matrix) + + def test_h_sparse_csr_format(self): + """h_sparse(format='csr') returns CSR matrix""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + h_sparse = graph.h_sparse(format='csr') + self.assertIsInstance(h_sparse, scipy.sparse.csr_matrix) + + def test_h_sparse_default_is_csc(self): + """h_sparse() default format is CSC""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + h_sparse = graph.h_sparse() + self.assertIsInstance(h_sparse, scipy.sparse.csc_matrix) + + def test_h_sparse_invalid_format_raises(self): + """h_sparse with invalid format raises ValueError""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + with self.assertRaises(Exception): + graph.h_sparse(format='invalid') + + def test_h_sparse_shape(self): + """h_sparse returns matrix with correct shape""" + graph = hg.HyperGraph() + + a, b, c = graph.new_variables([1, 2, 3]) + result = a * b + b * c + graph.compute(result) + + h_sparse = graph.h_sparse() + assert_equal(h_sparse.shape, (3, 3)) + + def test_h_sparse_multiplication(self): + """h_sparse values match dense h for multiplication""" + graph = hg.HyperGraph() + + a, b = graph.new_variables([5, 6]) + result = a * b + graph.compute(result) + + h_dense = graph.h(full=True) + h_sparse = graph.h_sparse(full=True) + + assert_array_almost_equal(h_sparse.toarray(), h_dense) + + def test_num_variables(self): + """num_variables property returns correct count""" + graph = hg.HyperGraph() + + assert_equal(graph.num_variables, 0) + + graph.new_variable(1) + assert_equal(graph.num_variables, 1) + + graph.new_variables([2, 3, 4]) + assert_equal(graph.num_variables, 4) + + if __name__ == '__main__': unittest.main()