feat: project named derivatives onto an explicit order

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.

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<index> positions(const std::vector<std::string> &names) const {
+    std::vector<index> 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<index>(it - m_d.begin())
+                           : -1);
+    }
+
+    return result;
+  }
+
   HYPERJET_INLINE static bool before(const std::pair<std::string, TScalar> &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<SScalar>
+  variables(const std::vector<std::pair<std::string, Scalar>> &values) {
+    std::vector<SScalar> 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<Scalar> g(const std::vector<std::string> &names) const {
+    const auto p = positions(names);
+
+    std::vector<Scalar> 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<Scalar> hm(const std::vector<std::string> &names) const
+    requires(TOrder == 2)
+  {
+    const auto p = positions(names);
+    const index n = length(names);
+    const index k = length(m_d);
+
+    std::vector<Scalar> result(static_cast<std::size_t>(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);

python/src/bind_sscalar.cpp

@@ -29,10 +29,55 @@ template <typename T> 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<std::string> &names) {
+            const auto g = self.g(names);
+
+            return py::array_t<typename T::Scalar>(
+                {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<std::pair<std::string, typename T::Scalar>> terms;
+        terms.reserve(values.size());
+
+        for (const auto &[name, value] : values) {
+          terms.emplace_back(name.template cast<std::string>(),
+                             value.template cast<typename T::Scalar>());
+        }
+
+        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<std::string> &names) {
+              const auto h = self.hm(names);
+              const auto n = hj::length(names);
+
+              py::array_t<typename T::Scalar> result({n, n});
+
+              std::copy(h.begin(), h.end(), result.mutable_data());
+
+              return result;
+            },
+            "names"_a);
   }
 
   // methods: arithmetic operations

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",

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))