From 40992b891f7dcb854430b06612705f74ffe7434b Mon Sep 17 00:00:00 2001 From: Thomas Date: Mon, 27 Jul 2026 21:17:06 +0200 Subject: [PATCH] feat: project named derivatives onto an explicit order Named derivatives were a dead end for anything that wants arrays. There was no way from an SScalar to a gradient vector or a Hessian matrix, and hj.d/hj.dd made it worse: they check for .g and .hm, a named scalar has neither, so they silently returned an empty array. A wrong answer rather than a missing feature. g(names) and hm(names) project onto an explicit, ordered list of names. The order comes from the caller, which is the point: values carrying different names project onto one shared layout without any of them being padded or remapped first. That is what an assembly step needs, and it is where named variables earn their keep over indexed ones -- no index bookkeeping between the local computation and the global system. A name the value does not carry reads as zero, the same rule d() and dd() already follow. hj.d and hj.dd take the same list via names=, and without it they now raise TypeError for named scalars instead of returning empty arrays. The indexed path is unchanged. variables(dict) builds a set in one call, in the order the dict gives, which is the order to hand back to g() and hm(). Lookups are batched: one binary search per requested name rather than one per matrix entry, so a k-by-k Hessian costs k log n + k^2 instead of k^2 log n. The mutation check on the new tests found one gap worth recording. The guard that skips unknown names in hm() can be removed without any assertion failing -- at(k, -1, j) reads before the buffer, and the garbage happened to be the zero the test expects. UBSan aborts on it, so the sanitizer job is what pins that guard, not the value comparisons. --- README.md | 31 +++++++ include/hyperjet/hyperjet.h | 80 ++++++++++++++++ python/src/bind_sscalar.cpp | 49 +++++++++- python/src/main.cpp | 57 ++++++++++-- python/tests/test_SScalar.py | 124 +++++++++++++++++++++++++ test/src/test_sscalar_second_order.cpp | 95 +++++++++++++++++++ 6 files changed, 426 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index ba236c3..e6e4997 100644 --- a/README.md +++ b/README.md @@ -165,6 +165,37 @@ The naming follows `DScalar`/`DDScalar`: one letter per order, so `SScalar` carr The Hessian is **dense over the names the value carries** — it is not sparse within that set. Any nonlinear function contributes the full outer product of its own gradient, so the triangle fills up after a couple of operations; what stays sparse is the set of names itself. It is stored as the upper triangle in the same packing `DDScalar` uses, and a first-order `SScalar` pays nothing for it: `sizeof` is unchanged at 32 bytes against 56 for `SSScalar`. +#### Getting back to arrays + +Named derivatives are only useful if they reach code that wants arrays. `g(names)` and `hm(names)` project onto an explicit, ordered list of names — the caller fixes the order, and a name the value does not carry reads as zero, just as `d` and `dd` do: + +```python +x, y, z = hj.SSScalar.variables({"x": 0.5, "y": 0.4, "z": 0.3}) + +f = x * y + +f.g(["x", "y", "z"]) # z is not in f, so it is zero +>>> array([0.4, 0.5, 0. ]) +f.hm(["x", "y", "z"]).shape +>>> (3, 3) +``` + +That the order comes from the caller is the point: values carrying *different* names project onto one shared layout, without any of them being padded or remapped first. This is what an assembly step needs, and it is where named variables beat indexed ones — no index bookkeeping between the local computation and the global system. + +The module-level `hj.d` and `hj.dd` take the same list: + +```python +values = [x * y, y * z] # the two carry different names + +hj.d(values, names=["x", "y", "z"]) +>>> array([[0.4, 0.5, 0. ], + [0. , 0.3, 0.4]]) +hj.dd(values, names=["x", "y", "z"]).shape +>>> (2, 3, 3) +``` + +Without `names`, `hj.d` and `hj.dd` raise `TypeError` for named scalars: a gradient does not exist until an order is chosen, and returning an empty array would be a wrong answer rather than a missing feature. + Two further differences from `DDScalar` worth knowing: - **No serialization.** `copy`, `deepcopy` and `pickle` raise `TypeError`; the `DDScalar` types support all three. diff --git a/include/hyperjet/hyperjet.h b/include/hyperjet/hyperjet.h index c1d1ac4..574a302 100644 --- a/include/hyperjet/hyperjet.h +++ b/include/hyperjet/hyperjet.h @@ -1800,6 +1800,24 @@ class SScalar { return (2 * k - 1 - i) * i / 2 + j; } + // Position of each requested name in the storage, or -1 where the value does + // not carry it. One lookup per name, so a projection costs k log n rather + // than a lookup per matrix entry. + std::vector positions(const std::vector &names) const { + std::vector result; + result.reserve(names.size()); + + for (const auto &name : names) { + const auto it = std::lower_bound(m_d.begin(), m_d.end(), name, before); + + result.push_back(it != m_d.end() && it->first == name + ? static_cast(it - m_d.begin()) + : -1); + } + + return result; + } + HYPERJET_INLINE static bool before(const std::pair &e, const std::string &name) { return e.first < name; @@ -1972,6 +1990,21 @@ class SScalar { return from_terms(value, Terms{{name, Scalar(1)}}); } + // Build a set of variables in one call. The order given here is the order + // g() and hm() project onto when handed the same names, so a caller can use + // one list for both. + static std::vector + variables(const std::vector> &values) { + std::vector result; + result.reserve(values.size()); + + for (const auto &[name, value] : values) { + result.push_back(variable(name, value)); + } + + return result; + } + private: // methods // Takes the storage as it is, already sorted. A second constructor would make // SScalar(f, {{"x", 1.0}}) ambiguous between Data and Terms. @@ -2014,6 +2047,26 @@ class SScalar { return r; } + // Projection onto an explicit, ordered list of names -- the bridge from + // named derivatives to code that wants arrays, an assembly step for + // instance. A name the value does not carry reads as zero, exactly as d() + // and dd() do, so a local value projects onto a larger global set without + // having to be padded first. The caller fixes the order, which is what makes + // the result composable: two values with different name sets project onto + // the same layout. + std::vector g(const std::vector &names) const { + const auto p = positions(names); + + std::vector result; + result.reserve(names.size()); + + for (const auto i : p) { + result.push_back(i < 0 ? Scalar(0) : m_d[i].second); + } + + return result; + } + Scalar dd(const std::string &a, const std::string &b) const requires(TOrder == 2) { @@ -2027,6 +2080,33 @@ class SScalar { return m_h[at(length(m_d), ia - m_d.begin(), ib - m_d.begin())]; } + // The Hessian over an explicit, ordered list of names, row-major and + // symmetric. Same zero-for-unknown rule as g(). + std::vector hm(const std::vector &names) const + requires(TOrder == 2) + { + const auto p = positions(names); + const index n = length(names); + const index k = length(m_d); + + std::vector result(static_cast(n * n), Scalar(0)); + + for (index i = 0; i < n; i++) { + if (p[i] < 0) { + continue; + } + + for (index j = 0; j < n; j++) { + if (p[j] < 0) { + continue; + } + + result[i * n + j] = m_h[at(k, p[i], p[j])]; + } + } + + return result; + } Scalar eval(const Data &d) const { Scalar r(m_f); diff --git a/python/src/bind_sscalar.cpp b/python/src/bind_sscalar.cpp index bf23342..d91e082 100644 --- a/python/src/bind_sscalar.cpp +++ b/python/src/bind_sscalar.cpp @@ -29,10 +29,55 @@ template void bind_sscalar(pybind11::module &m, const char *name) { .def("abs", &T::abs) .def("d", &T::d, "variable"_a) .def("eval", &T::eval, "d"_a) - .def("names", &T::names); + .def("names", &T::names) + .def( + // Projection onto an explicit, ordered list of names. Returns an + // array so the result drops straight into array-consuming code; a + // name the value does not carry reads as zero. + "g", + [](const T &self, const std::vector &names) { + const auto g = self.g(names); + + return py::array_t( + {hj::length(g)}, {sizeof(typename T::Scalar)}, g.data()); + }, + "names"_a); + + // static methods + cls.def_static( + "variables", + [](const py::dict &values) { + // A dict keeps insertion order, and that order is what g() and hm() + // project onto when handed the same names. + std::vector> terms; + terms.reserve(values.size()); + + for (const auto &[name, value] : values) { + terms.emplace_back(name.template cast(), + value.template cast()); + } + + return T::variables(terms); + }, + "values"_a); if constexpr (T::order() == 2) { - cls.def("dd", &T::dd, "a"_a, "b"_a); + cls.def("dd", &T::dd, "a"_a, "b"_a) + .def( + // The Hessian over an explicit, ordered list of names, as an + // n-by-n array. + "hm", + [](const T &self, const std::vector &names) { + const auto h = self.hm(names); + const auto n = hj::length(names); + + py::array_t result({n, n}); + + std::copy(h.begin(), h.end(), result.mutable_data()); + + return result; + }, + "names"_a); } // methods: arithmetic operations diff --git a/python/src/main.cpp b/python/src/main.cpp index 6dc85a3..b8536d3 100644 --- a/python/src/main.cpp +++ b/python/src/main.cpp @@ -37,14 +37,55 @@ PYBIND11_MODULE(hyperjet, m) { auto global = py::dict(); global["np"] = numpy; - m.attr("f") = py::eval( - "np.vectorize(lambda v: v.f if hasattr(v, 'f') else v)", global); - m.attr("d") = py::eval("np.vectorize(lambda v: v.g if hasattr(v, 'g') else " - "np.zeros((0)), signature='()->(n)')", - global); - m.attr("dd") = py::eval("np.vectorize(lambda v: v.hm() if hasattr(v, 'hm') " - "else np.zeros((0, 0)), signature='()->(n,m)')", - global); + // Indexed scalars carry their own order, so d() and dd() can read the + // gradient and Hessian directly. Named scalars have no canonical order -- + // each value carries whichever names it picked up -- so the caller has to + // supply one via names=, which is also what makes a set of values with + // differing names project onto a single shared layout. + py::exec(R"PY( +def _indexed_g(v): + if hasattr(v, "names"): + raise TypeError( + "hj.d needs an order for named scalars: pass names=[...]" + ) + return v.g if hasattr(v, "g") else np.zeros(0) + + +def _indexed_hm(v): + if hasattr(v, "names"): + raise TypeError( + "hj.dd needs an order for named scalars: pass names=[...]" + ) + return v.hm() if hasattr(v, "hm") else np.zeros((0, 0)) + + +_f = np.vectorize(lambda v: v.f if hasattr(v, "f") else v) +_g = np.vectorize(_indexed_g, signature="()->(n)") +_hm = np.vectorize(_indexed_hm, signature="()->(n,m)") +_named_g = np.vectorize(lambda v, n: v.g(n), signature="()->(n)", excluded={1}) +_named_hm = np.vectorize(lambda v, n: v.hm(n), signature="()->(n,m)", excluded={1}) + + +def f(values): + return _f(values) + + +def d(values, names=None): + if names is None: + return _g(values) + return _named_g(values, list(names)) + + +def dd(values, names=None): + if names is None: + return _hm(values) + return _named_hm(values, list(names)) +)PY", + global); + + m.attr("f") = global["f"]; + m.attr("d") = global["d"]; + m.attr("dd") = global["dd"]; m.def( "variables", diff --git a/python/tests/test_SScalar.py b/python/tests/test_SScalar.py index 0b54e87..f971226 100644 --- a/python/tests/test_SScalar.py +++ b/python/tests/test_SScalar.py @@ -345,3 +345,127 @@ def test_hessian_of_an_unknown_name_is_zero(): assert_allclose(r.dd("x", "z"), 0) assert_allclose(r.dd("z", "z"), 0) + + +# projection onto an explicit list of names +# +# Named derivatives are only useful if they reach code that wants arrays. The +# caller supplies the order, which is what lets values carrying different names +# project onto one shared layout -- an assembly step, for instance. + + +def test_variables_builds_a_set_in_dict_order(): + x, y, z = hj.SSScalar.variables({"x": 0.5, "y": 0.4, "z": 0.3}) + + assert_allclose([x.f, y.f, z.f], [0.5, 0.4, 0.3]) + + # each carries only its own name + assert_equal([x.names(), y.names(), z.names()], [["x"], ["y"], ["z"]]) + assert_allclose(x.d("x"), 1) + assert_allclose(x.d("y"), 0) + + +def test_g_uses_the_given_order(): + x, y, z = hj.SSScalar.variables({"x": 0.5, "y": 0.4, "z": 0.3}) + + r = x * y + z + + assert_allclose(r.g(["x", "y", "z"]), [r.d("x"), r.d("y"), r.d("z")]) + + # the caller's order, not the storage order + assert_allclose(r.g(["z", "y", "x"]), [r.d("z"), r.d("y"), r.d("x")]) + + # a name the value does not carry reads as zero, as with d + assert_allclose(r.g(["w", "x"]), [0.0, r.d("x")]) + + assert_equal(r.g([]).shape, (0,)) + + +def test_hm_is_square_and_symmetric(): + x, y, z = hj.SSScalar.variables({"x": 0.5, "y": 0.4, "z": 0.3}) + names = ["x", "y", "z"] + + h = (x * y + z * z).sqrt().hm(names) + + assert_equal(h.shape, (3, 3)) + assert_allclose(h, h.T) + + r = (x * y + z * z).sqrt() + + for i, a in enumerate(names): + for j, b in enumerate(names): + assert_allclose(h[i, j], r.dd(a, b)) + + +def test_hm_zeroes_unknown_names(): + x, y = hj.SSScalar.variables({"x": 0.5, "y": 0.4}) + + h = (x * y).hm(["x", "w", "y"]) + + assert_allclose(h[1, :], 0) + assert_allclose(h[:, 1], 0) + assert_allclose(h[0, 2], (x * y).dd("x", "y")) + + +def test_projection_returns_a_copy(): + x, y = hj.SSScalar.variables({"x": 0.5, "y": 0.4}) + + r = x * y + g = r.g(["x", "y"]) + + g[0] = 100.0 + + # writing to the result does not touch the scalar + assert_allclose(r.d("x"), 0.4) + assert_allclose(r.g(["x", "y"]), [0.4, 0.5]) + + +def test_first_order_has_no_hm(): + assert hasattr(hj.SScalar, "g") + assert not hasattr(hj.SScalar, "hm") + assert hasattr(hj.SSScalar, "hm") + + +def test_utilities_need_an_order_for_named_scalars(): + x, y = hj.SSScalar.variables({"x": 0.5, "y": 0.4}) + + values = [x * y, x + y] + + # the value is unambiguous, so f needs nothing + assert_allclose(hj.f(values), [0.2, 0.9]) + + # a gradient does not exist without an order -- silently returning an empty + # array would be the wrong answer, not a missing feature + with pytest.raises(TypeError, match="names"): + hj.d(values) + + with pytest.raises(TypeError, match="names"): + hj.dd(values) + + +def test_utilities_project_a_set_onto_one_layout(): + x, y, z = hj.SSScalar.variables({"x": 0.5, "y": 0.4, "z": 0.3}) + names = ["x", "y", "z"] + + # the two values carry different names + values = [x * y, y * z] + + g = hj.d(values, names=names) + + assert_equal(g.shape, (2, 3)) + assert_allclose(g[0], [0.4, 0.5, 0.0]) + assert_allclose(g[1], [0.0, 0.3, 0.4]) + + h = hj.dd(values, names=names) + + assert_equal(h.shape, (2, 3, 3)) + assert_allclose(h[0], (x * y).hm(names)) + assert_allclose(h[1], (y * z).hm(names)) + + +def test_utilities_still_work_for_indexed_scalars(): + values = [v**2 for v in hj.variables([1.0, 2.0, 3.0])] + + assert_allclose(hj.f(values), [1.0, 4.0, 9.0]) + assert_equal(hj.d(values).shape, (3, 3)) + assert_equal(hj.dd(values).shape, (3, 3, 3)) diff --git a/test/src/test_sscalar_second_order.cpp b/test/src/test_sscalar_second_order.cpp index ce984d6..2e051d2 100644 --- a/test/src/test_sscalar_second_order.cpp +++ b/test/src/test_sscalar_second_order.cpp @@ -281,4 +281,99 @@ TEST_CASE("SScalar second order: a variable has no curvature") { CHECK(u().dd("x", "q") == doctest::Approx(0.0)); } +// Projection onto an explicit list of names -- the bridge to code that wants +// arrays. The values come from d() and dd(), pinned symbolically above; what +// is new here is the layout and the zero-for-unknown rule. + +TEST_CASE("SScalar g projects onto the given order") { + const auto r = vx() * vy() + vz(); + + // the caller's order, not the storage order + const auto g = r.g({"z", "y", "x"}); + + CHECK(hyperjet::length(g) == 3); + CHECK(g[0] == doctest::Approx(r.d("z"))); + CHECK(g[1] == doctest::Approx(r.d("y"))); + CHECK(g[2] == doctest::Approx(r.d("x"))); + + // a name the value does not carry reads as zero + const auto padded = r.g({"w", "x"}); + + CHECK(padded[0] == doctest::Approx(0.0)); + CHECK(padded[1] == doctest::Approx(r.d("x"))); + + // an empty list is empty, not an error + CHECK(hyperjet::length(r.g({})) == 0); +} + +TEST_CASE("SScalar hm is row-major and symmetric") { + using std::sqrt; + + const auto r = sqrt(vx() * vy() + vz() * vz()); + const std::vector names{"x", "y", "z"}; + + const auto h = r.hm(names); + + CHECK(hyperjet::length(h) == 9); + + for (hyperjet::index i = 0; i < 3; i++) { + for (hyperjet::index j = 0; j < 3; j++) { + CHECK(h[i * 3 + j] == doctest::Approx(r.dd(names[i], names[j]))); + CHECK(h[i * 3 + j] == doctest::Approx(h[j * 3 + i])); + } + } +} + +TEST_CASE("SScalar hm zeroes rows and columns of unknown names") { + const auto r = vx() * vy(); + const std::vector names{"x", "w", "y"}; + + const auto h = r.hm(names); + + const auto entry = [&](const hyperjet::index i, const hyperjet::index j) { + return h[i * 3 + j]; + }; + + // w is unknown, so its row and column are zero throughout + for (hyperjet::index i = 0; i < 3; i++) { + CHECK(entry(1, i) == doctest::Approx(0.0)); + CHECK(entry(i, 1) == doctest::Approx(0.0)); + } + + // the remaining entries keep the positions the caller chose + CHECK(entry(0, 2) == doctest::Approx(r.dd("x", "y"))); + CHECK(entry(2, 0) == doctest::Approx(r.dd("x", "y"))); +} + +TEST_CASE("SScalar projects a local value onto a larger set") { + // Two values carrying different names project onto one shared layout, which + // is what an assembly step needs. + const std::vector global{"x", "y", "z"}; + + const auto a = vx() * vy(); + const auto b = vy() * vz(); + + CHECK(a.g(global)[2] == doctest::Approx(0.0)); // a does not know z + CHECK(b.g(global)[0] == doctest::Approx(0.0)); // b does not know x + + CHECK(hyperjet::length(a.hm(global)) == 9); + CHECK(hyperjet::length(b.hm(global)) == 9); +} + +TEST_CASE("SScalar variables builds a set in the given order") { + const auto v = S::variables({{"a", 1.0}, {"b", 2.0}, {"c", 3.0}}); + + CHECK(hyperjet::length(v) == 3); + + CHECK(v[0].f() == doctest::Approx(1.0)); + CHECK(v[1].f() == doctest::Approx(2.0)); + CHECK(v[2].f() == doctest::Approx(3.0)); + + CHECK(v[0].d("a") == doctest::Approx(1.0)); + CHECK(v[0].d("b") == doctest::Approx(0.0)); + + // each carries only its own name + CHECK(v[1].names() == std::vector{"b"}); +} + } // namespace