From 9c09e19ae895d27448e53a4d418f31b1ed2ed664 Mon Sep 17 00:00:00 2001 From: Thomas Date: Mon, 27 Jul 2026 17:44:23 +0200 Subject: [PATCH] refactor: store the SScalar derivatives sorted instead of in a hash map A multiply of two single-entry values cost six allocations, eval copied the name on every iteration, and binary and ternary copied the whole map. The storage becomes a vector of name and derivative sorted by name, so combining two values is a linear merge rather than a sequence of hash lookups and node allocations. Measured on an expression chain, -O3: names before after 2 416 ns, 12 allocs 86 ns, 2 allocs 4.8x 4 2029 ns, 53 allocs 280 ns, 6 allocs 7.2x 12 14271 ns, 375 allocs 1546 ns, 22 allocs 9.2x This replaces the design proposed for this step. Interning the names behind a shared table was going to be the approach, and the allocation accounting did not survive scrutiny: with different name sets the merge and the remap of both operands cost no less than the map, and in (x*y)/(x-y) the tables are equal in content but not in identity, so the fast path would never have fired. Both designs were prototyped and measured before choosing. The sorted vector is not a detour for second order either. The position of a name in the vector is its index, which is what a dense triangular Hessian would be laid out over, and the merge already produces the union in sorted order. Data stays the map. It is the interface type for construction and eval, which is what callers and in particular a Python dict naturally provide; the vector is internal. A second constructor taking the internal type would have made SScalar(f, {{"x", 1.0}}) ambiguous between the two, so the operations go through a named factory instead. Three simplifications fall out. The four in-place operators delegate to their binary form, which puts the merge in one place and makes aliased operands correct without the guards added in #76. The derivatives are now iterated in a defined order, so operator<< and repr are deterministic instead of unspecified -- the tests and the README are tightened accordingly. And clang-tidy drops from 50 findings to 41, because the braceless single statements it complained about were in the code that went away. Also replaces the int coefficients in operator+ and operator- with Scalar, which matters for a float scalar type. The public API is unchanged: 35 bound methods as before, SScalar(f=, d={...}), size, d and eval identical. That the characterization from #84 carries the rewrite without a single expected value being touched is the evidence that doing it first was right. C++ 112 test cases, 1421 assertions, Debug and Release with -Werror and under clang with ASan+UBSan Python 248 passed clang-format and ruff format clean --- README.md | 2 +- include/hyperjet/hyperjet.h | 225 +++++++++++++++++------------------ python/tests/test_SScalar.py | 9 +- test/src/test_sscalar.cpp | 11 +- 4 files changed, 113 insertions(+), 134 deletions(-) diff --git a/README.md b/README.md index 56fd79c..6c37f94 100644 --- a/README.md +++ b/README.md @@ -143,7 +143,7 @@ len(u * v) # the product knows both names Two further differences from `DDScalar` worth knowing: - **No serialization.** `copy`, `deepcopy` and `pickle` raise `TypeError`; the `DDScalar` types support all three. -- **The order of terms in `repr` is unspecified**, because the derivatives live in an unordered map. Only the value comes first. +- **The gradient is stored sorted by name**, so `repr` and iteration are deterministic, and combining two values is a linear merge rather than a sequence of hash lookups. `eval(d)` contracts the gradient with a displacement given per name. Names missing from the displacement contribute nothing, and names the scalar does not know are ignored: diff --git a/include/hyperjet/hyperjet.h b/include/hyperjet/hyperjet.h index 7872737..8bc89b0 100644 --- a/include/hyperjet/hyperjet.h +++ b/include/hyperjet/hyperjet.h @@ -12,6 +12,7 @@ #include // string #include // conditional #include // unordered_map +#include // move, pair #include // vector namespace hyperjet { @@ -1744,55 +1745,79 @@ template class SScalar { public: // types using Type = SScalar; using Scalar = TScalar; + + // The interface type for construction and eval. A map is what callers, and + // in particular a Python dict, naturally provide. using Data = std::unordered_map; +private: // types + // The storage. Sorted by name, so a combination of two values is a linear + // merge rather than a sequence of hash lookups, and the position of a name is + // its index -- which is what a second-order Hessian would be laid out over. + using Terms = std::vector>; + private: // variables // The default constructor is defaulted, so without an initializer m_f would // be indeterminate and reading it undefined. Scalar m_f{}; - Data m_d; + Terms m_d; private: // methods + HYPERJET_INLINE static bool before(const std::pair &e, + const std::string &name) { + return e.first < name; + } + template - HYPERJET_INLINE static Data unary(const Data &a, const TDa &da) { - Data r{a}; + HYPERJET_INLINE static Terms unary(const Terms &a, const TDa &da) { + Terms r{a}; + // scaling the derivatives leaves the order intact for (auto &d : r) d.second *= da; return r; } + // da * a + db * b over the union of the names, in one pass template - HYPERJET_INLINE static Data binary(const Data &a, const Data &b, - const TDa &da, const TDb &db) { - Data r{a}; + HYPERJET_INLINE static Terms binary(const Terms &a, const Terms &b, + const TDa &da, const TDb &db) { + Terms r; + r.reserve(a.size() + b.size()); + + std::size_t i = 0; + std::size_t j = 0; + + while (i < a.size() && j < b.size()) { + if (a[i].first < b[j].first) { + r.emplace_back(a[i].first, da * a[i].second); + i++; + } else if (b[j].first < a[i].first) { + r.emplace_back(b[j].first, db * b[j].second); + j++; + } else { + r.emplace_back(a[i].first, da * a[i].second + db * b[j].second); + i++; + j++; + } + } - for (auto &d : r) - d.second *= da; + for (; i < a.size(); i++) + r.emplace_back(a[i].first, da * a[i].second); - for (auto &d : b) - r[d.first] += db * d.second; + for (; j < b.size(); j++) + r.emplace_back(b[j].first, db * b[j].second); return r; } template - HYPERJET_INLINE static Data ternary(const Data &a, const Data &b, - const Data &c, const TDa &da, - const TDb &db, const TDc &dc) { - Data r{a}; - - for (auto &d : r) - d.second *= da; - - for (auto &d : b) - r[d.first] += db * d.second; - - for (auto &d : c) - r[d.first] += dc * d.second; - - return r; + HYPERJET_INLINE static Terms ternary(const Terms &a, const Terms &b, + const Terms &c, const TDa &da, + const TDb &db, const TDc &dc) { + // two merges rather than a three-way one; only hypot of three uses this + return binary(binary(a, b, da, db), c, Scalar(1), dc); } public: // constructors @@ -1800,23 +1825,37 @@ template class SScalar { SScalar(const Scalar f) : m_f(f), m_d() {} - SScalar(const Scalar f, const Data &d) : m_f(f), m_d(d) {} + SScalar(const Scalar f, const Data &d) : m_f(f), m_d(d.begin(), d.end()) { + std::sort(m_d.begin(), m_d.end(), + [](const auto &a, const auto &b) { return a.first < b.first; }); + } static SScalar constant(const Scalar value) { return SScalar(value, Data()); } static SScalar variable(const std::string &name, const Scalar value) { - return SScalar(value, {{name, Scalar(1)}}); + return from_terms(value, Terms{{name, Scalar(1)}}); + } + +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. + HYPERJET_INLINE static SScalar from_terms(const Scalar f, Terms d) { + SScalar r; + r.m_f = f; + r.m_d = std::move(d); + return r; } +public: public: // methods index size() const { return m_d.size(); } Scalar f() const { return m_f; } Scalar d(const std::string &variable) const { - const auto it = m_d.find(variable); + const auto it = std::lower_bound(m_d.begin(), m_d.end(), variable, before); - if (it == m_d.end()) + if (it == m_d.end() || it->first != variable) return Scalar(0); return it->second; @@ -1875,10 +1914,10 @@ template class SScalar { Type operator+(const Type &b) const { const auto f = m_f + b.f(); - const auto da = 1; - const auto db = 1; + const auto da = Scalar(1); + const auto db = Scalar(1); - return Type(f, binary(m_d, b.m_d, da, db)); + return from_terms(f, binary(m_d, b.m_d, da, db)); } Type operator+(const Scalar b) const { @@ -1891,14 +1930,9 @@ template class SScalar { friend Type operator+(const Scalar a, const Type &b) { return b + a; } - Type &operator+=(const Type &b) { - m_f += b.m_f; - - for (const auto &d : b.m_d) - m_d[d.first] += d.second; - - return *this; - } + // The in-place forms delegate, which keeps the merge in one place and makes + // aliased operands correct without a guard of their own. + Type &operator+=(const Type &b) { return *this = *this + b; } Type &operator+=(const Scalar &b) { m_f += b; @@ -1910,10 +1944,10 @@ template class SScalar { Type operator-(const Type &b) const { const auto f = m_f - b.f(); - const auto da = 1; - const auto db = -1; + const auto da = Scalar(1); + const auto db = Scalar(-1); - return Type(f, binary(m_d, b.m_d, da, db)); + return from_terms(f, binary(m_d, b.m_d, da, db)); } Type operator-(const Scalar b) const { return -b + *this; } @@ -1926,14 +1960,7 @@ template class SScalar { return result; } - Type &operator-=(const Type &b) { - m_f -= b.m_f; - - for (const auto &d : b.m_d) - m_d[d.first] -= d.second; - - return *this; - } + Type &operator-=(const Type &b) { return *this = *this - b; } Type &operator-=(const Scalar &b) { m_f -= b; @@ -1948,7 +1975,7 @@ template class SScalar { const auto da = b.f(); const auto db = m_f; - return Type(f, binary(m_d, b.m_d, da, db)); + return from_terms(f, binary(m_d, b.m_d, da, db)); } Type operator*(const Scalar b) const { @@ -1964,25 +1991,7 @@ template class SScalar { friend Type operator*(const Scalar a, const Type &b) { return b * a; } - Type &operator*=(const Type &b) { - // aliased operands would be read again after being overwritten below - if (this == &b) { - return *this = *this * b; - } - - const Scalar da = b.m_f; - const Scalar db = m_f; - - m_f *= b.m_f; - - for (auto &d : m_d) - d.second *= da; - - for (auto &d : b.m_d) - m_d[d.first] += db * d.second; - - return *this; - } + Type &operator*=(const Type &b) { return *this = *this * b; } Type &operator*=(const Scalar &b) { m_f *= b; @@ -2004,7 +2013,7 @@ template class SScalar { const auto da = tmp; const auto db = -m_f / pow(b.f(), 2); - return Type(f, binary(m_d, b.m_d, da, db)); + return from_terms(f, binary(m_d, b.m_d, da, db)); } Type operator/(const Scalar b) const { return Scalar(1) / b * (*this); } @@ -2015,30 +2024,10 @@ template class SScalar { const auto f = a / b.f(); const auto db = -a / pow(b.f(), 2); - return Type(f, unary(b.m_d, db)); + return from_terms(f, unary(b.m_d, db)); } - Type &operator/=(const Type &b) { - using std::pow; - - // aliased operands would be read again after being overwritten below - if (this == &b) { - return *this = *this / b; - } - - const auto da = 1 / b.m_f; - const auto db = -m_f / pow(b.m_f, 2); - - m_f /= b.m_f; - - for (auto &d : m_d) - d.second *= da; - - for (auto &d : b.m_d) - m_d[d.first] += db * d.second; - - return *this; - } + Type &operator/=(const Type &b) { return *this = *this / b; } Type &operator/=(const Scalar &b) { operator*=(1 / b); @@ -2054,7 +2043,7 @@ template class SScalar { const auto f = pow(m_f, b); const auto da = b * pow(m_f, b - Scalar(1)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type sqrt() const { @@ -2063,7 +2052,7 @@ template class SScalar { const auto f = sqrt(m_f); const auto da = Scalar(1) / (Scalar(2) * f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type cbrt() const { @@ -2072,14 +2061,14 @@ template class SScalar { const auto f = cbrt(m_f); const auto da = Scalar(1) / (Scalar(3) * f * f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type reciprocal() const { const auto f = Scalar(1) / m_f; const auto da = -f * f; - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } // trigonometric @@ -2091,7 +2080,7 @@ template class SScalar { const auto f = cos(m_f); const auto da = -sin(m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type sin() const { @@ -2101,7 +2090,7 @@ template class SScalar { const auto f = sin(m_f); const auto da = cos(m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type tan() const { @@ -2110,7 +2099,7 @@ template class SScalar { const auto f = tan(m_f); const auto da = f * f + Scalar(1); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type acos() const { @@ -2122,7 +2111,7 @@ template class SScalar { const auto f = acos(m_f); const auto da = -Scalar(1) / sqrt(tmp); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type asin() const { @@ -2134,7 +2123,7 @@ template class SScalar { const auto f = asin(m_f); const auto da = Scalar(1) / sqrt(tmp); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type atan() const { @@ -2143,7 +2132,7 @@ template class SScalar { const auto f = atan(m_f); const auto da = Scalar(1) / (m_f * m_f + Scalar(1)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type atan2(const Type &b) const { @@ -2155,7 +2144,7 @@ template class SScalar { const auto da = b.m_f / tmp; const auto db = -m_f / tmp; - return Type(f, binary(m_d, b.m_d, da, db)); + return from_terms(f, binary(m_d, b.m_d, da, db)); } static Type hypot(const Type &a, const Type &b) { @@ -2165,7 +2154,7 @@ template class SScalar { const auto da = a.m_f / f; const auto db = b.m_f / f; - return Type(f, binary(a.m_d, b.m_d, da, db)); + return from_terms(f, binary(a.m_d, b.m_d, da, db)); } static Type hypot(const Type &a, const Type &b, const Type &c) { @@ -2177,7 +2166,7 @@ template class SScalar { const auto db = b.m_f / f; const auto dc = c.m_f / f; - return Type(f, ternary(a.m_d, b.m_d, c.m_d, da, db, dc)); + return from_terms(f, ternary(a.m_d, b.m_d, c.m_d, da, db, dc)); } // hyperbolic @@ -2189,7 +2178,7 @@ template class SScalar { const auto f = cosh(m_f); const auto da = sinh(m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type sinh() const { @@ -2199,7 +2188,7 @@ template class SScalar { const auto f = sinh(m_f); const auto da = cosh(m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type tanh() const { @@ -2208,7 +2197,7 @@ template class SScalar { const auto f = tanh(m_f); const auto da = 1 - f * f; - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type acosh() const { @@ -2218,7 +2207,7 @@ template class SScalar { const auto f = acosh(m_f); const auto da = 1 / (sqrt(m_f - 1) * sqrt(m_f + 1)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type asinh() const { @@ -2228,7 +2217,7 @@ template class SScalar { const auto f = asinh(m_f); const auto da = 1 / sqrt(1 + m_f * m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type atanh() const { @@ -2238,7 +2227,7 @@ template class SScalar { const auto f = atanh(m_f); const auto da = 1 / (1 - m_f * m_f); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } // exponents and logarithms @@ -2249,7 +2238,7 @@ template class SScalar { const auto f = exp(m_f); const auto da = f; - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type log() const { @@ -2258,7 +2247,7 @@ template class SScalar { const auto f = log(m_f); const auto da = 1 / m_f; - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type log(const TScalar base) const { @@ -2267,7 +2256,7 @@ template class SScalar { const auto f = log(m_f) / log(base); const auto da = 1 / (m_f * log(base)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type log2() const { @@ -2277,7 +2266,7 @@ template class SScalar { const auto f = log2(m_f); const auto da = 1 / (m_f * log(2)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } Type log10() const { @@ -2287,7 +2276,7 @@ template class SScalar { const auto f = log10(m_f); const auto da = 1 / (m_f * log(10)); - return Type(f, unary(m_d, da)); + return from_terms(f, unary(m_d, da)); } // abs diff --git a/python/tests/test_SScalar.py b/python/tests/test_SScalar.py index b28ffd8..cab8ca3 100644 --- a/python/tests/test_SScalar.py +++ b/python/tests/test_SScalar.py @@ -101,13 +101,8 @@ def test_properties(): def test_repr(): - # The derivatives live in an unordered map, so the order of the terms is - # unspecified. Only the value comes first and every term appears once. - text = repr(hj.SScalar(f=3, d={"x": 1, "y": 6})) - - assert text.startswith("3") - assert "+1*dx" in text - assert "+6*dy" in text + # The derivatives are stored sorted by name, so the order is deterministic. + assert_equal(repr(hj.SScalar(f=3, d={"y": 6, "x": 1})), "3 +1*dx +6*dy") # negative derivatives are printed without a plus sign assert_equal(repr(hj.SScalar(f=1, d={"x": -2})), "1 -2*dx") diff --git a/test/src/test_sscalar.cpp b/test/src/test_sscalar.cpp index 0e6e420..17feeb2 100644 --- a/test/src/test_sscalar.cpp +++ b/test/src/test_sscalar.cpp @@ -185,17 +185,12 @@ TEST_CASE("SScalar size counts the stored derivatives") { } TEST_CASE("SScalar printing") { - // The derivatives live in an unordered_map, so the order of the terms is - // unspecified. Only the value comes first, and every term appears once. + // The derivatives are stored sorted by name, so the order of the terms is + // deterministic. std::stringstream out; out << s2; - const std::string text = out.str(); - CHECK(text.starts_with("4")); - CHECK(text.find("*dx") != std::string::npos); - CHECK(text.find("*dy") != std::string::npos); - CHECK(text.find("+7") != std::string::npos); - CHECK(text.find("+1") != std::string::npos); + CHECK(out.str() == "4 +7*dx +1*dy"); // negative derivatives are printed without a plus sign std::stringstream negative;