diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 94a7615..c8b31bc 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -20,8 +20,13 @@ jobs: - uses: astral-sh/setup-uv@v6 + # The file list is discovered rather than spelled out: a hand-maintained + # list silently stops covering a new directory, and an entry whose glob + # matches nothing makes clang-format fail on the literal pattern. The .in + # templates are excluded on purpose -- they hold CMake @SUBSTITUTIONS@ and + # are not valid C++ until configured. - name: Check C++ formatting - run: uv run --only-group dev clang-format --dry-run --Werror include/hyperjet/*.h python/src/*.h python/src/*.cpp test/src/*.cpp benchmark/src/*.cpp + run: uv run --only-group dev clang-format --dry-run --Werror $(git ls-files '*.h' '*.hpp' '*.cpp') - name: Check Python formatting run: uv run --only-group dev ruff format --check python/tests/ diff --git a/README.md b/README.md index e6e4997..edf8ffc 100644 --- a/README.md +++ b/README.md @@ -94,12 +94,12 @@ HyperJet provides two families of scalar types: Stores derivatives in dense arrays. Supports first-order (gradient only) and second-order (gradient + Hessian) derivatives. Available in **static** variants with a compile-time-fixed number of variables, and a **dynamic** variant for arbitrary sizes. -| Python type | Order | Variables | C++ type | -|-------------|-------|-----------|----------| -| `DScalar` | 1 | dynamic | `DDScalar<1, double>` | -| `DDScalar` | 2 | dynamic | `DDScalar<2, double>` | -| `D3Scalar` | 1 | 3 (static) | `DDScalar<1, double, 3>` | -| `DD3Scalar` | 2 | 3 (static) | `DDScalar<2, double, 3>` | +| Python type | Order | Variables | C++ type | +| ----------- | ----- | ---------- | ------------------------ | +| `DScalar` | 1 | dynamic | `DDScalar<1, double>` | +| `DDScalar` | 2 | dynamic | `DDScalar<2, double>` | +| `D3Scalar` | 1 | 3 (static) | `DDScalar<1, double, 3>` | +| `DD3Scalar` | 2 | 3 (static) | `DDScalar<2, double, 3>` | Static variants (`D0Scalar`–`D16Scalar`, `DD0Scalar`–`DD16Scalar`) avoid heap allocation and enable better compiler optimization. The dynamic variants (`DScalar`, `DDScalar`) accept any number of variables at runtime. @@ -180,7 +180,7 @@ 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. +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: @@ -214,6 +214,77 @@ u.eval({"w": 100.0}) # w is unknown and ignored >>> 3.0 ``` +## Performance + +Derivatives of one function, computed every way a caller might reasonably choose. The function is the strain energy of a cable with unit-spaced nodes and transverse displacements `x_i`, so each term is the squared elongation of one segment: + +``` +E(x) = Σ_i ( sqrt(1 + (x_{i+1} - x_i)²) - 1 )² +``` + +Median nanoseconds per evaluation, seven repetitions. The spread within a run is under 1 %, and the AD libraries reproduce to about 1 % across independent runs and across a fresh build. Two kinds of figure are looser: those in the single digits sit at the timer's resolution, and the finite-difference rows vary by up to 5 % between process starts. + +### Second order — value, gradient and dense Hessian + +| | n=3 | n=6 | n=12 | n=24 | +| ------------------------------ | ------ | ------- | ------- | ---------- | +| **HyperJet, static size** | **15** | **101** | **921** | **12_214** | +| Eigen `AutoDiffScalar`, nested | 18 | 328 | 2_888 | 22_105 | +| autodiff `dual2nd` | 127 | 375 | 2_425 | 17_679 | +| HyperJet, dynamic size | 644 | 1928 | 6_296 | 22_777 | +| central finite differences | 513 | 2659 | 9_968 | 55_604 | + +HyperJet is fastest at every size, by **3.3× at n=6 and 2.6× at n=12** against whichever of the other two libraries is faster there, narrowing to 1.5× at n=24 and 1.2× at n=3. + +Finite differences are 4.6× to 34× slower _and_ approximate — the worst deviation from the exact Hessian rises to 8e-8 at n=24. HyperJet's own dynamic variant allocates per operation and costs 43× the static one at n=3, falling to 1.9× at n=24 where the arithmetic dominates the allocation. + +### First order — value and gradient + +Ceres joins here. `ceres::Jet` is first order only, which is not a limitation but what a solver needs: it wants Jacobians. Nesting `Jet` inside `Jet` to force it into the table above would benchmark a configuration no Ceres user writes. + +| | n=3 | n=6 | n=12 | n=24 | +| -------------------------- | --- | ---- | ---- | ------- | +| HyperJet, static size | 4 | 12 | 59 | **244** | +| Eigen `AutoDiffScalar` | 3 | 10 | 57 | 245 | +| `ceres::Jet` | 3 | 11 | 81 | 304 | +| autodiff `dual` | 76 | 94 | 213 | 725 | +| central finite differences | 173 | 397 | 819 | 1946 | +| HyperJet, dynamic size | 664 | 1614 | 3465 | 6889 | + +**Here there is no winner.** HyperJet and `AutoDiffScalar` are level from n=12 upward — 244 against 245 at n=24 — and both are ahead of `Jet`, which costs 1.25× as much at n=24. Below n=12 the three sit within a few nanoseconds of each other, which is the timer's resolution rather than a result. + +Gradients alone are not where a hyper-dual library can distinguish itself: every contender stores a value and one dense derivative vector and does the same arithmetic on it. The second-order table is where the designs diverge. + +One caveat about both tables, and it is not specific to HyperJet. Every one of these libraries carries its derivatives inside the scalar and returns a whole scalar by value from every operation — 200 bytes at n=24, first order — so a chained expression depends on the compiler keeping those temporaries in registers rather than spilling them. Whether it does is decided by an inlining threshold on the _enclosing_ expression, and falling on the wrong side of it costs all of them about the same: + +| first order, n=24 | inlined | not inlined | penalty | +| ---------------------- | ------- | ----------- | ------- | +| HyperJet | 246 | 383 | 1.56× | +| Eigen `AutoDiffScalar` | 247 | 418 | 1.69× | +| `ceres::Jet` | 307 | 500 | 1.63× | + +The ordering survives either way. The benchmark forces the inlining uniformly, because otherwise the threshold rather than the library decides the comparison — the compiler happened to inline for some contenders and not others at the same call site. + +In a hot loop it is worth writing `acc += term` rather than `acc = acc + term`: the compound operators work in place, with no returned temporary at all. On the second-order n=24 case that alone is worth 8 %. + +### What follows from this + +- **Second order is where HyperJet earns its place.** That is what it was built for: one pass produces the full Hessian, where `Jet` and `AutoDiffScalar` need nesting and `dual2nd` re-evaluates the function n(n+1)/2 times. +- **If you only need gradients, take whichever library you already have** — Eigen is level, Ceres is close, and Ceres brings a solver with it. +- **Take the static types when the variable count is known at compile time.** The dynamic variant trades most of the performance away. + +The benchmark verifies that every contender produces the same value, gradient and Hessian before it times anything — a fast contender computing the wrong thing would otherwise read as a win. All the AD libraries agree bit-for-bit; only finite differences deviate. + +Reproduce with: + +``` +cmake -Sbenchmark -Bbuild/benchmark -DCMAKE_BUILD_TYPE=Release +cmake --build build/benchmark --target Compare +./build/benchmark/Compare --benchmark_repetitions=7 --benchmark_report_aggregates_only=true +``` + +Measured on an Apple M1 Pro, macOS 26.5, Apple clang 21.0.0, `-O3 -DNDEBUG`, against Eigen 3.4.0, autodiff v1.1.2 and Ceres 2.2.0. Absolute numbers will differ on other machines; the ratios are the point. + ## Validation Arguments coming from Python are validated. Out-of-range Hessian indices raise `IndexError`; inconsistent sizes raise `RuntimeError` — `eval` with the wrong number of values, `set_hm` with a mismatched shape, a negative size, a variable index outside the gradient, or padding below the current size. @@ -273,14 +344,14 @@ a.nbytes # 3 x 80 bytes, contiguous Arithmetic and the mathematical functions then run as compiled loops instead of a Python object loop. Measured over 10 000 elements of `DD3Scalar`: -| | `dtype=object` | dtype | | -|---|---|---|---| -| `a + b` | 349.5 ns | 2.5 ns | 138× | -| `a * b` | 354.3 ns | 3.9 ns | 92× | -| `np.sqrt(a)` | 349.2 ns | 2.8 ns | 123× | -| `np.sin(a)` | 366.4 ns | 7.0 ns | 52× | -| `np.sum(a)` | 353.1 ns | 3.8 ns | 94× | -| `a @ b` | 673.0 ns | 2.7 ns | 246× | +| | `dtype=object` | dtype | | +| ------------ | -------------- | ------ | ---- | +| `a + b` | 349.5 ns | 2.5 ns | 138× | +| `a * b` | 354.3 ns | 3.9 ns | 92× | +| `np.sqrt(a)` | 349.2 ns | 2.8 ns | 123× | +| `np.sin(a)` | 366.4 ns | 7.0 ns | 52× | +| `np.sum(a)` | 353.1 ns | 3.8 ns | 94× | +| `a @ b` | 673.0 ns | 2.7 ns | 246× | The **dynamic** variants (`DScalar`, `DDScalar`) hold a `std::vector`, so they have no fixed element size and stay object arrays. @@ -331,13 +402,13 @@ The library requires a C++23-capable compiler (GCC ≥ 14, Clang ≥ 18, MSVC All of the scalar types support: -| Category | Functions | -|----------|-----------| -| Arithmetic | `+`, `-`, `*`, `/`, `pow`, `abs`, `reciprocal` | +| Category | Functions | +| ------------- | ---------------------------------------------------- | +| Arithmetic | `+`, `-`, `*`, `/`, `pow`, `abs`, `reciprocal` | | Trigonometric | `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `atan2` | -| Hyperbolic | `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh` | -| Exponential | `exp`, `log`, `log2`, `log10` | -| Other | `sqrt`, `cbrt`, `hypot` | +| Hyperbolic | `sinh`, `cosh`, `tanh`, `asinh`, `acosh`, `atanh` | +| Exponential | `exp`, `log`, `log2`, `log10` | +| Other | `sqrt`, `cbrt`, `hypot` | ## Utility Functions diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 9158e8b..0a29e3f 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -32,6 +32,37 @@ endif() CPMAddPackage("gh:/pybind/pybind11@2.13.6") +# For the comparison benchmark only: another forward-mode AD library to measure +# against. Header-only, so DOWNLOAD_ONLY is enough. +CPMAddPackage( + NAME autodiff + GITHUB_REPOSITORY autodiff/autodiff + GIT_TAG v1.1.2 + DOWNLOAD_ONLY YES +) + +# Ceres for the first-order comparison. ceres::Jet is first order only, and +# jet.h needs nothing but its own include tree and Eigen -- no glog, no Abseil, +# so the whole solver does not have to be built. +CPMAddPackage( + NAME ceres + GITHUB_REPOSITORY ceres-solver/ceres-solver + GIT_TAG 2.2.0 + DOWNLOAD_ONLY YES +) + +if(ceres_ADDED) + add_library(ceres_jet INTERFACE IMPORTED) + target_include_directories(ceres_jet INTERFACE ${ceres_SOURCE_DIR}/include) +endif() + +if(autodiff_ADDED) + add_library(autodiff INTERFACE IMPORTED) + target_include_directories(autodiff INTERFACE ${autodiff_SOURCE_DIR}) + # autodiff gates its Eigen overloads on this rather than detecting Eigen + target_compile_definitions(autodiff INTERFACE AUTODIFF_EIGEN_FOUND) +endif() + if(TEST_INSTALLED_VERSION) find_package(HyperJet REQUIRED) else() @@ -51,3 +82,18 @@ target_compile_definitions(Benchmark PRIVATE ) set_target_properties(Benchmark PROPERTIES CXX_STANDARD 23) + +# --- Comparison benchmark +# +# A separate target: it pulls in autodiff, which the micro-benchmarks have no +# use for, and it has its own main() that verifies the contenders agree before +# timing anything. + +add_executable(Compare ${CMAKE_CURRENT_SOURCE_DIR}/compare/compare.cpp) + +target_link_libraries(Compare PRIVATE benchmark::benchmark Eigen autodiff + ceres_jet hyperjet::hyperjet) + +target_compile_definitions(Compare PRIVATE -DHYPERJET_EXCEPTIONS) + +set_target_properties(Compare PROPERTIES CXX_STANDARD 23) diff --git a/benchmark/compare/compare.cpp b/benchmark/compare/compare.cpp new file mode 100644 index 0000000..c48c87e --- /dev/null +++ b/benchmark/compare/compare.cpp @@ -0,0 +1,760 @@ +// Derivatives of one non-trivial function, computed every way a caller might +// reasonably choose -- first order, then first and second order together. +// +// The micro-benchmarks in ../src measure single operations against nothing. +// This measures the thing a caller actually wants -- value, gradient and full +// Hessian of a real function -- against the alternatives someone choosing a +// library would consider. +// +// The function is the strain energy of a cable with unit-spaced nodes and +// transverse displacements x_i: +// +// E(x) = sum_i ( sqrt(1 + (x_{i+1} - x_i)^2) - 1 )^2 +// +// Each term is the squared elongation of one segment. It is nonlinear, its +// Hessian is not constant, and it is defined for any n, which is what makes it +// usable across all four contenders and all sizes. +// +// The Hessian is banded, but every contender here computes all n^2 entries +// regardless -- dense forward mode does not exploit structure -- so the shape +// does not favour any of them. +// +// Note that the two orders are separate comparisons. ceres::Jet is first order +// only, which is not a limitation of Ceres but what a solver needs: it wants +// Jacobians. Nesting Jet inside Jet to force it into the second-order table +// would benchmark a configuration no Ceres user writes. + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include // cbrt, sqrt, abs +#include // printf +#include // abort +#include // string +#include +#include + +// The compiler decides per instantiation whether to inline term() below, and it +// decided differently for different contenders: HyperJet's body is larger, so +// clang left it out of line and every call returned a whole scalar -- 200 bytes +// at n=24 -- through the stack. That cost 1.5x and had nothing to do with the +// libraries, only with where an inlining threshold happened to fall. Forcing it +// uniformly measures the arithmetic instead of the threshold. +#if defined _MSC_VER +#define COMPARE_INLINE __forceinline +#else +#define COMPARE_INLINE __attribute__((always_inline)) +#endif + +namespace hj = hyperjet; + +namespace { + +// The function, written once. Each contender supplies an accessor, so all four +// differentiate the same expression rather than four transcriptions of it. +template auto energy(const int n, TAt at) { + using std::sqrt; + + // The accessor's plain scalar type. Every intermediate is spelled out as this + // type on purpose: both Eigen's AutoDiffScalar and autodiff's dual are + // expression-template types, so `auto d = at(i + 1) - at(i)` would keep an + // expression that refers to operands which are gone by the time it is read, + // and returning `auto` from term() would do the same. HyperJet has no + // expression templates and does not care either way. + // + // Starting the accumulator from a literal zero is no help: not every + // contender's scalar is constructible from a double. + using Scalar = std::decay_t; + + const auto term = [&](const int i) COMPARE_INLINE -> Scalar { + const Scalar d = at(i + 1) - at(i); + const Scalar s = sqrt(1.0 + d * d) - 1.0; + + return s * s; + }; + + Scalar e = term(0); + + for (int i = 1; i + 1 < n; i++) { + e = e + term(i); + } + + return e; +} + +std::vector point(const int n) { + std::vector x(static_cast(n)); + + for (int i = 0; i < n; i++) { + // an arbitrary but non-degenerate configuration; a straight cable would + // make every segment's elongation zero and flatter the nonlinear terms + x[static_cast(i)] = 0.1 * std::sin(0.7 * (i + 1)); + } + + return x; +} + +// The result all four have to agree on: value, gradient, dense Hessian. +// +// The caller allocates it once and every contender fills it in place. Filling a +// fresh one inside the timing loop would put three vector allocations into +// every measurement -- at n=3 that is the same order as the whole computation, +// so it would measure std::vector rather than the differentiation. +struct Result { + double f{}; + std::vector g; + std::vector h; // row-major, n by n + + void resize(const int n) { + g.resize(static_cast(n)); + h.resize(static_cast(n * n)); + } + + // For the first-order contenders. deviation() walks whatever is there, so + // leaving h empty simply takes it out of the comparison. + void resize_gradient(const int n) { + g.resize(static_cast(n)); + h.clear(); + } +}; + +// --- HyperJet, static size + +template +void hyperjet_static(const std::vector &x, Result &r) { + using S = hj::DDScalar<2, double, TSize>; + + std::array values{}; + + for (int i = 0; i < TSize; i++) { + values[static_cast(i)] = x[static_cast(i)]; + } + + const auto v = S::variables(values); + const auto e = energy(TSize, [&](const int i) { return v[i]; }); + + r.f = e.f(); + + for (int i = 0; i < TSize; i++) { + r.g[static_cast(i)] = e.g(i); + + for (int j = 0; j < TSize; j++) { + r.h[static_cast(i * TSize + j)] = e.h(i, j); + } + } +} + +// --- HyperJet, dynamic size + +void hyperjet_dynamic(const std::vector &x, Result &r) { + using S = hj::DDScalar<2, double>; + + const auto n = static_cast(x.size()); + const auto v = S::variables(x); + const auto e = energy(n, [&](const int i) { return v[i]; }); + + r.f = e.f(); + + for (int i = 0; i < n; i++) { + r.g[static_cast(i)] = e.g(i); + + for (int j = 0; j < n; j++) { + r.h[static_cast(i * n + j)] = e.h(i, j); + } + } +} + +// --- Eigen AutoDiffScalar, nested for second order +// +// Eigen ships first-order forward mode. Second order comes from nesting: the +// outer scalar's derivatives are themselves first-order AD scalars, so the +// outer gradient carries the Hessian rows. Seeding that takes O(n^2) writes +// before the function is even called. + +template +void eigen_nested(const std::vector &x, Result &r) { + using Inner = Eigen::AutoDiffScalar>; + using Outer = Eigen::AutoDiffScalar>; + + Eigen::Matrix v; + + for (int i = 0; i < TSize; i++) { + v(i).value().value() = x[static_cast(i)]; + v(i).value().derivatives() = Eigen::Matrix::Unit(i); + + for (int j = 0; j < TSize; j++) { + v(i).derivatives()(j).value() = i == j ? 1.0 : 0.0; + v(i).derivatives()(j).derivatives().setZero(); + } + } + + const auto e = energy(TSize, [&](const int i) { return v(i); }); + + r.f = e.value().value(); + + for (int i = 0; i < TSize; i++) { + r.g[static_cast(i)] = e.value().derivatives()(i); + + for (int j = 0; j < TSize; j++) { + r.h[static_cast(i * TSize + j)] = + e.derivatives()(i).derivatives()(j); + } + } +} + +// --- autodiff, dual2nd +// +// autodiff's hessian() seeds one pair (i, j) at a time and re-evaluates the +// function for each, so it calls the function n(n+1)/2 times rather than once. +// That is the library's design, not a misuse of it: dual2nd carries a single +// second-order direction. + +void autodiff_dual2nd(const std::vector &x, Result &r) { + const auto n = static_cast(x.size()); + + autodiff::VectorXdual2nd v(n); + + for (int i = 0; i < n; i++) { + v(i) = x[static_cast(i)]; + } + + const auto fn = [n](const auto &values) { + return energy(n, [&](const int i) { return values(i); }); + }; + + autodiff::dual2nd u; + Eigen::VectorXd g; + const Eigen::MatrixXd h = + autodiff::hessian(fn, autodiff::wrt(v), autodiff::at(v), u, g); + + r.f = static_cast(u); + + for (int i = 0; i < n; i++) { + r.g[static_cast(i)] = g(i); + + for (int j = 0; j < n; j++) { + r.h[static_cast(i * n + j)] = h(i, j); + } + } +} + +// --- Central finite differences +// +// The baseline. No library, no types -- and no exact answer either, which is +// the point of including it: the deviation is reported alongside the timing. + +void finite_differences(const std::vector &x, Result &r) { + const auto n = static_cast(x.size()); + + // cbrt(eps) is the usual choice for second derivatives: it balances + // truncation against cancellation + const double h = std::cbrt(std::numeric_limits::epsilon()); + + const auto eval = [](const std::vector &values) { + const auto m = static_cast(values.size()); + + return energy( + m, [&](const int i) { return values[static_cast(i)]; }); + }; + + const auto shifted = [&](const int i, const double di, const int j, + const double dj) { + auto values = x; + values[static_cast(i)] += di; + values[static_cast(j)] += dj; + + return eval(values); + }; + + r.f = eval(x); + + for (int i = 0; i < n; i++) { + const auto fp = shifted(i, h, i, 0.0); + const auto fm = shifted(i, -h, i, 0.0); + + r.g[static_cast(i)] = (fp - fm) / (2.0 * h); + r.h[static_cast(i * n + i)] = (fp - 2.0 * r.f + fm) / (h * h); + + for (int j = i + 1; j < n; j++) { + const auto pp = shifted(i, h, j, h); + const auto pm = shifted(i, h, j, -h); + const auto mp = shifted(i, -h, j, h); + const auto mm = shifted(i, -h, j, -h); + + const auto value = (pp - pm - mp + mm) / (4.0 * h * h); + + r.h[static_cast(i * n + j)] = value; + r.h[static_cast(j * n + i)] = value; + } + } +} + +// === First order: value and gradient only === +// +// A separate comparison, because this is the regime Ceres and Eigen were built +// for and where most callers actually live. + +// --- HyperJet, static size + +template +void hyperjet_static_g(const std::vector &x, Result &r) { + using S = hj::DDScalar<1, double, TSize>; + + std::array values{}; + + for (int i = 0; i < TSize; i++) { + values[static_cast(i)] = x[static_cast(i)]; + } + + const auto v = S::variables(values); + const auto e = energy(TSize, [&](const int i) { return v[i]; }); + + r.f = e.f(); + + for (int i = 0; i < TSize; i++) { + r.g[static_cast(i)] = e.g(i); + } +} + +// --- HyperJet, dynamic size + +void hyperjet_dynamic_g(const std::vector &x, Result &r) { + using S = hj::DDScalar<1, double>; + + const auto n = static_cast(x.size()); + const auto v = S::variables(x); + const auto e = energy(n, [&](const int i) { return v[i]; }); + + r.f = e.f(); + + for (int i = 0; i < n; i++) { + r.g[static_cast(i)] = e.g(i); + } +} + +// --- Ceres Jet +// +// First order by design. The gradient lives in a fixed-size Eigen vector, so +// there is no dynamic counterpart to measure. + +template void ceres_jet(const std::vector &x, Result &r) { + using J = ceres::Jet; + + std::array v; + + for (int i = 0; i < TSize; i++) { + auto &jet = v[static_cast(i)]; + + jet.a = x[static_cast(i)]; + jet.v.setZero(); + jet.v[i] = 1.0; + } + + const auto e = energy(TSize, [&](const int i) { return v[i]; }); + + r.f = e.a; + + for (int i = 0; i < TSize; i++) { + r.g[static_cast(i)] = e.v[i]; + } +} + +// --- Eigen AutoDiffScalar, single level + +template +void eigen_first_order(const std::vector &x, Result &r) { + using AD = Eigen::AutoDiffScalar>; + + Eigen::Matrix v; + + for (int i = 0; i < TSize; i++) { + v(i).value() = x[static_cast(i)]; + v(i).derivatives() = Eigen::Matrix::Unit(i); + } + + const auto e = energy(TSize, [&](const int i) { return v(i); }); + + r.f = e.value(); + + for (int i = 0; i < TSize; i++) { + r.g[static_cast(i)] = e.derivatives()(i); + } +} + +// --- autodiff, dual + +void autodiff_dual(const std::vector &x, Result &r) { + const auto n = static_cast(x.size()); + + autodiff::VectorXdual v(n); + + for (int i = 0; i < n; i++) { + v(i) = x[static_cast(i)]; + } + + const auto fn = [n](const auto &values) { + return energy(n, [&](const int i) { return values(i); }); + }; + + autodiff::dual u; + const Eigen::VectorXd g = + autodiff::gradient(fn, autodiff::wrt(v), autodiff::at(v), u); + + r.f = static_cast(u); + + for (int i = 0; i < n; i++) { + r.g[static_cast(i)] = g(i); + } +} + +// --- Central finite differences, gradient only + +void finite_differences_g(const std::vector &x, Result &r) { + const auto n = static_cast(x.size()); + + // cbrt(eps) suits second derivatives; a gradient alone does better with + // sqrt(eps), and giving finite differences its best shot is the fair way to + // include it + const double h = std::sqrt(std::numeric_limits::epsilon()); + + const auto eval = [](const std::vector &values) { + const auto m = static_cast(values.size()); + + return energy( + m, [&](const int i) { return values[static_cast(i)]; }); + }; + + r.f = eval(x); + + for (int i = 0; i < n; i++) { + auto plus = x; + auto minus = x; + + plus[static_cast(i)] += h; + minus[static_cast(i)] -= h; + + r.g[static_cast(i)] = (eval(plus) - eval(minus)) / (2.0 * h); + } +} + +// --- Verification +// +// Timing a contender that computes the wrong thing is worthless, so every +// contender is checked against HyperJet before any measurement runs. + +double deviation(const Result &a, const Result &b) { + auto worst = std::abs(a.f - b.f); + + for (std::size_t i = 0; i < a.g.size(); i++) { + worst = std::max(worst, std::abs(a.g[i] - b.g[i])); + } + + for (std::size_t i = 0; i < a.h.size(); i++) { + worst = std::max(worst, std::abs(a.h[i] - b.h[i])); + } + + return worst; +} + +template void verify() { + const auto x = point(TSize); + + Result reference; + reference.resize(TSize); + hyperjet_static(x, reference); + + const auto run = [&](auto fn) { + Result r; + r.resize(TSize); + fn(x, r); + + return deviation(reference, r); + }; + + const struct { + const char *name; + double worst; + double tolerance; + } checks[] = { + {"HyperJet dynamic", run(hyperjet_dynamic), 0.0}, + {"Eigen nested", run(eigen_nested), 1e-15}, + {"autodiff dual2nd", run(autodiff_dual2nd), 1e-15}, + // finite differences are the odd one out: no tolerance this tight can + // hold, and the number itself is the interesting result + {"finite differences", run(finite_differences), 1e-4}, + }; + + for (const auto &check : checks) { + std::printf(" n=%-3d %-20s worst deviation %.3e %s\n", TSize, check.name, + check.worst, + check.worst <= check.tolerance ? "" : "<-- FAILED"); + + if (check.worst > check.tolerance) { + std::printf("\nA contender disagrees with HyperJet; timings would be " + "meaningless.\n"); + std::abort(); + } + } +} + +template void verify_first_order() { + const auto x = point(TSize); + + Result reference; + reference.resize_gradient(TSize); + hyperjet_static_g(x, reference); + + const auto run = [&](auto fn) { + Result r; + r.resize_gradient(TSize); + fn(x, r); + + return deviation(reference, r); + }; + + const struct { + const char *name; + double worst; + double tolerance; + } checks[] = { + {"HyperJet dynamic", run(hyperjet_dynamic_g), 0.0}, + {"Ceres Jet", run(ceres_jet), 1e-16}, + {"Eigen AutoDiffScalar", run(eigen_first_order), 1e-16}, + {"autodiff dual", run(autodiff_dual), 1e-16}, + {"finite differences", run(finite_differences_g), 1e-7}, + }; + + for (const auto &check : checks) { + std::printf(" n=%-3d %-20s worst deviation %.3e %s\n", TSize, check.name, + check.worst, + check.worst <= check.tolerance ? "" : "<-- FAILED"); + + if (check.worst > check.tolerance) { + std::printf("\nA contender disagrees with HyperJet; timings would be " + "meaningless.\n"); + std::abort(); + } + } +} + +// --- Benchmarks + +template void hyperjet_static_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize(TSize); + + for (auto _ : state) { + // The input is loop-invariant and the computation a pure function of it, + // so without this the whole thing can be hoisted out of the loop and the + // measurement becomes fiction. + benchmark::DoNotOptimize(x); + hyperjet_static(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void hyperjet_dynamic_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize(TSize); + + for (auto _ : state) { + // The input is loop-invariant and the computation a pure function of it, + // so without this the whole thing can be hoisted out of the loop and the + // measurement becomes fiction. + benchmark::DoNotOptimize(x); + hyperjet_dynamic(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void eigen_nested_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize(TSize); + + for (auto _ : state) { + // The input is loop-invariant and the computation a pure function of it, + // so without this the whole thing can be hoisted out of the loop and the + // measurement becomes fiction. + benchmark::DoNotOptimize(x); + eigen_nested(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void autodiff_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize(TSize); + + for (auto _ : state) { + // The input is loop-invariant and the computation a pure function of it, + // so without this the whole thing can be hoisted out of the loop and the + // measurement becomes fiction. + benchmark::DoNotOptimize(x); + autodiff_dual2nd(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void finite_differences_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize(TSize); + + for (auto _ : state) { + // The input is loop-invariant and the computation a pure function of it, + // so without this the whole thing can be hoisted out of the loop and the + // measurement becomes fiction. + benchmark::DoNotOptimize(x); + finite_differences(x, r); + benchmark::DoNotOptimize(r); + } +} + +// first order + +template void hyperjet_static_g_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + hyperjet_static_g(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void hyperjet_dynamic_g_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + hyperjet_dynamic_g(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void ceres_jet_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + ceres_jet(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void eigen_first_order_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + eigen_first_order(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void autodiff_dual_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + autodiff_dual(x, r); + benchmark::DoNotOptimize(r); + } +} + +template void finite_differences_g_bm(benchmark::State &state) { + auto x = point(TSize); + + Result r; + r.resize_gradient(TSize); + + for (auto _ : state) { + benchmark::DoNotOptimize(x); + finite_differences_g(x, r); + benchmark::DoNotOptimize(r); + } +} + +} // namespace + +#define COMPARE_FIRST_ORDER(size) \ + BENCHMARK_TEMPLATE(hyperjet_static_g_bm, size); \ + BENCHMARK_TEMPLATE(hyperjet_dynamic_g_bm, size); \ + BENCHMARK_TEMPLATE(ceres_jet_bm, size); \ + BENCHMARK_TEMPLATE(eigen_first_order_bm, size); \ + BENCHMARK_TEMPLATE(autodiff_dual_bm, size); \ + BENCHMARK_TEMPLATE(finite_differences_g_bm, size); + +COMPARE_FIRST_ORDER(3) +COMPARE_FIRST_ORDER(6) +COMPARE_FIRST_ORDER(12) +COMPARE_FIRST_ORDER(24) + +#define COMPARE(size) \ + BENCHMARK_TEMPLATE(hyperjet_static_bm, size); \ + BENCHMARK_TEMPLATE(hyperjet_dynamic_bm, size); \ + BENCHMARK_TEMPLATE(eigen_nested_bm, size); \ + BENCHMARK_TEMPLATE(autodiff_bm, size); \ + BENCHMARK_TEMPLATE(finite_differences_bm, size); + +COMPARE(3) +COMPARE(6) +COMPARE(12) +COMPARE(24) + +int main(int argc, char **argv) { + std::printf("Verifying first order:\n"); + + verify_first_order<3>(); + verify_first_order<6>(); + verify_first_order<12>(); + verify_first_order<24>(); + + std::printf("\nVerifying first and second order:\n"); + + verify<3>(); + verify<6>(); + verify<12>(); + verify<24>(); + + std::printf("\n"); + + benchmark::Initialize(&argc, argv); + benchmark::RunSpecifiedBenchmarks(); + benchmark::Shutdown(); + + return 0; +} diff --git a/include/hyperjet/hyperjet.h b/include/hyperjet/hyperjet.h index 574a302..b8ab16a 100644 --- a/include/hyperjet/hyperjet.h +++ b/include/hyperjet/hyperjet.h @@ -164,6 +164,16 @@ class DDScalar { static HYPERJET_INLINE Scalar operator()(const Scalar b) { return b; } }; + // Filling the Hessian triangle in one pass halves its memory traffic, but it + // also shortens the gradient pass to the gradient alone, and below a handful + // of variables that one long vectorised pass is worth more than the traffic + // saved. Measured crossover at second order, Apple M1 Pro, clang 21: n=6 is + // 14% slower fused, n=8 3% faster, n=16 35% faster. + // + // size() is a compile-time constant for the static variants, so the branch + // folds away there. + static constexpr index hessian_fusion_threshold = 8; + template static constexpr bool is_tag_v = std::is_same_v || std::is_same_v || @@ -179,68 +189,87 @@ class DDScalar { } } - template + // The chain rule term da * H_a has the same shape as the gradient term, so + // one pass over the whole data array fills the Hessian slots as well. Only a + // curvature term makes the triangle need a loop of its own. + // + // r never aliases a: the in-place operators copy their operand first, which + // is what lets the triangle read a's gradient after r's has been written. + template HYPERJET_INLINE void unary(const Data &a, const Scalar f, const TDa da, const TDaa daa, Data &r) const noexcept { - const index n = data_length(); - r[0] = f; - if constexpr (TOrder < 1 || std::is_same()) { + if constexpr (TOrder < 1 || std::is_same_v) { return; - } + } else if constexpr (TOrder < 2 || std::is_same_v) { + const index n = data_length(); - for (index i = 1; i < n; i++) { - if constexpr (TIncrement) { - r[i] += mul(da, a[i]); - } else { + for (index i = 1; i < n; i++) { r[i] = mul(da, a[i]); } - } + } else if (size() < hessian_fusion_threshold) { + const index n = data_length(); - if constexpr (TOrder < 2 || std::is_same()) { - return; - } + for (index i = 1; i < n; i++) { + r[i] = mul(da, a[i]); + } - index k = 1 + size(); + index k = 1 + size(); - for (index i = 0; i < size(); i++) { - const auto ca = mul(daa, a[1 + i]); + for (index i = 0; i < size(); i++) { + const auto ca = mul(daa, a[1 + i]); - for (index j = i; j < size(); j++) { - r[k++] += ca * a[1 + j]; + for (index j = i; j < size(); j++) { + r[k++] += ca * a[1 + j]; + } + } + } else { + for (index i = 1; i <= size(); i++) { + r[i] = mul(da, a[i]); + } + + index k = 1 + size(); + + for (index i = 0; i < size(); i++) { + const auto ca = mul(daa, a[1 + i]); + + for (index j = i; j < size(); j++) { + r[k] = mul(da, a[k]) + ca * a[1 + j]; + k++; + } } } } - template + // See unary() for why the Hessian is one pass and why aliasing is not a + // concern here. + template HYPERJET_INLINE void binary(const Data &a, const Data &b, const Scalar f, const TDa da, const TDb db, const TDaa daa, const TDab dab, const TDbb dbb, Data &r) const noexcept { - const index n = data_length(); - r[0] = f; if constexpr (TOrder < 1 || (std::is_same_v && std::is_same_v)) { return; - } else { + } else if constexpr (TOrder < 2 || (std::is_same_v && + std::is_same_v && + std::is_same_v)) { + const index n = data_length(); + for (index i = 1; i < n; i++) { - if constexpr (TIncrement) { - r[i] += mul(da, a[i]) + mul(db, b[i]); - } else { - r[i] = mul(da, a[i]) + mul(db, b[i]); - } + r[i] = mul(da, a[i]) + mul(db, b[i]); + } + } else if (size() < hessian_fusion_threshold) { + const index n = data_length(); + + for (index i = 1; i < n; i++) { + r[i] = mul(da, a[i]) + mul(db, b[i]); } - } - if constexpr (TOrder < 2 || - (std::is_same_v && std::is_same_v && - std::is_same_v)) { - return; - } else { index k = 1 + size(); for (index i = 0; i < size(); i++) { @@ -251,41 +280,59 @@ class DDScalar { r[k++] += ca * a[1 + j] + cb * b[1 + j]; } } + } else { + for (index i = 1; i <= size(); i++) { + r[i] = mul(da, a[i]) + mul(db, b[i]); + } + + index k = 1 + size(); + + for (index i = 0; i < size(); i++) { + const auto ca = mul(daa, a[1 + i]) + mul(dab, b[1 + i]); + const auto cb = mul(dab, a[1 + i]) + mul(dbb, b[1 + i]); + + for (index j = i; j < size(); j++) { + r[k] = mul(da, a[k]) + mul(db, b[k]) + ca * a[1 + j] + cb * b[1 + j]; + k++; + } + } } } - template + // See unary() for why the Hessian is one pass and why aliasing is not a + // concern here. + template HYPERJET_INLINE void ternary(const Data &a, const Data &b, const Data &c, const Scalar f, const TDa da, const TDb db, const TDc dc, const TDaa daa, const TDab dab, const TDac dac, const TDbb dbb, const TDbc dbc, const TDcc dcc, Data &r) const noexcept { - const index n = data_length(); - r[0] = f; if constexpr (TOrder < 1 || (std::is_same_v && std::is_same_v && std::is_same_v)) { return; - } else { + } else if constexpr (TOrder < 2 || (std::is_same_v && + std::is_same_v && + std::is_same_v && + std::is_same_v && + std::is_same_v && + std::is_same_v)) { + const index n = data_length(); + for (index i = 1; i < n; i++) { - if constexpr (TIncrement) { - r[i] += mul(da, a[i]) + mul(db, b[i]) + mul(dc, c[i]); - } else { - r[i] = mul(da, a[i]) + mul(db, b[i]) + mul(dc, c[i]); - } + r[i] = mul(da, a[i]) + mul(db, b[i]) + mul(dc, c[i]); + } + } else if (size() < hessian_fusion_threshold) { + const index n = data_length(); + + for (index i = 1; i < n; i++) { + r[i] = mul(da, a[i]) + mul(db, b[i]) + mul(dc, c[i]); } - } - if constexpr (TOrder < 2 || - (std::is_same_v && std::is_same_v && - std::is_same_v && std::is_same_v && - std::is_same_v && std::is_same_v)) { - return; - } else { index k = 1 + size(); for (index i = 0; i < size(); i++) { @@ -300,6 +347,27 @@ class DDScalar { r[k++] += ca * a[1 + j] + cb * b[1 + j] + cc * c[1 + j]; } } + } else { + for (index i = 1; i <= size(); i++) { + r[i] = mul(da, a[i]) + mul(db, b[i]) + mul(dc, c[i]); + } + + index k = 1 + size(); + + for (index i = 0; i < size(); i++) { + const auto ca = + mul(daa, a[1 + i]) + mul(dab, b[1 + i]) + mul(dac, c[1 + i]); + const auto cb = + mul(dab, a[1 + i]) + mul(dbb, b[1 + i]) + mul(dbc, c[1 + i]); + const auto cc = + mul(dac, a[1 + i]) + mul(dbc, b[1 + i]) + mul(dcc, c[1 + i]); + + for (index j = i; j < size(); j++) { + r[k] = mul(da, a[k]) + mul(db, b[k]) + mul(dc, c[k]) + ca * a[1 + j] + + cb * b[1 + j] + cc * c[1 + j]; + k++; + } + } } } @@ -334,6 +402,16 @@ class DDScalar { static_assert(!is_dynamic()); } + // Takes ownership instead of copying. The factories build their data as a + // temporary, and for the dynamic variant copying it in meant a second + // allocation and a second pass over the whole array. + DDScalar(Data &&data, const index size) + : m_size(size), m_data(std::move(data)) { + static_assert(0 < order() && order() <= 2); + + static_assert(is_dynamic()); + } + DDScalar(const Data &data, const index size) : m_size(size), m_data(data) { static_assert(0 < order() && order() <= 2); @@ -483,40 +561,38 @@ class DDScalar { } } - static Type empty() { + HYPERJET_INLINE static Type empty() { if constexpr (is_dynamic()) { - Data data(1); - Type result(data, 0); - return result; + // Size zero, not the default constructor's size one. + return Type(Data(1), 0); } else { - Data data; - Type result(data); - return result; + // The default constructor already leaves m_data indeterminate, which is + // what empty() means. Building a local Data and copying it in was both a + // wasted pass over the whole array -- 2.6 kB at order 2 with 24 + // variables -- and a read of indeterminate doubles, which is undefined. + return Type(); } } static Type empty(const index size) { if constexpr (is_dynamic()) { check_valid_size(size); - const index n = data_length_from_size(size); - const Data data(n); - Type result(data, size); - return result; + return Type(Data(data_length_from_size(size)), size); } else { check_valid_size(size); return empty(); } } - static Type zero() { + HYPERJET_INLINE static Type zero() { if constexpr (is_dynamic()) { - Data data(1); - Type result(data, 0); - return result; + return Type(Data(1), 0); } else { - Data data; - data.fill(0); - Type result(data); + // Filling in place rather than filling a local and copying it in: the + // copy was a second pass over the array, and variables() pays it once per + // variable. + Type result; + result.m_data.fill(0); return result; } } @@ -524,9 +600,7 @@ class DDScalar { static Type zero(const index size) { if constexpr (is_dynamic()) { check_valid_size(size); - const Data data(data_length_from_size(size), 0); - Type result(data, size); - return result; + return Type(Data(data_length_from_size(size), 0), size); } else { check_valid_size(size); return zero(); @@ -819,7 +893,7 @@ class DDScalar { // --- neg - Type operator-() const { + HYPERJET_INLINE Type operator-() const { Type result = Type::empty(size()); for (index i = 0; i < data_length(); i++) { @@ -831,7 +905,7 @@ class DDScalar { // --- add - Type operator+(const Type &b) const { + HYPERJET_INLINE Type operator+(const Type &b) const { check_equal_size(size(), b.size()); Type result = Type::empty(size()); @@ -845,12 +919,12 @@ class DDScalar { const auto dab = Zero(); const auto dbb = Zero(); - binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); + binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } - Type operator+(const Scalar b) const { + HYPERJET_INLINE Type operator+(const Scalar b) const { Type result = *this; result.f() += b; @@ -860,7 +934,7 @@ class DDScalar { friend Type operator+(const Scalar a, const Type &b) { return b + a; } - Type &operator+=(const Type &b) { + HYPERJET_INLINE Type &operator+=(const Type &b) { check_equal_size(size(), b.size()); for (index i = 0; i < length(m_data); i++) { @@ -870,7 +944,7 @@ class DDScalar { return *this; } - Type &operator+=(const Scalar &b) { + HYPERJET_INLINE Type &operator+=(const Scalar &b) { f() += b; return *this; @@ -878,7 +952,7 @@ class DDScalar { // --- sub - Type operator-(const Type &b) const { + HYPERJET_INLINE Type operator-(const Type &b) const { check_equal_size(size(), b.size()); Type result = Type::empty(size()); @@ -892,12 +966,12 @@ class DDScalar { const auto dab = Zero(); const auto dbb = Zero(); - binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); + binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } - Type operator-(const Scalar b) const { return -b + *this; } + HYPERJET_INLINE Type operator-(const Scalar b) const { return -b + *this; } friend Type operator-(const Scalar a, const Type &b) { Type result = Type::empty(b.size()); @@ -911,7 +985,7 @@ class DDScalar { return result; } - Type &operator-=(const Type &b) { + HYPERJET_INLINE Type &operator-=(const Type &b) { check_equal_size(size(), b.size()); for (index i = 0; i < length(m_data); i++) { @@ -921,7 +995,7 @@ class DDScalar { return *this; } - Type &operator-=(const Scalar &b) { + HYPERJET_INLINE Type &operator-=(const Scalar &b) { f() -= b; return *this; @@ -929,7 +1003,7 @@ class DDScalar { // --- mul - Type operator*(const Type &b) const { + HYPERJET_INLINE Type operator*(const Type &b) const { check_equal_size(size(), b.size()); Type result = Type::empty(size()); @@ -943,12 +1017,12 @@ class DDScalar { const auto dab = One(); const auto dbb = Zero(); - binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); + binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } - Type operator*(const Scalar b) const { + HYPERJET_INLINE Type operator*(const Scalar b) const { Type result = Type::empty(size()); for (index i = 0; i < data_length(); i++) { @@ -960,7 +1034,7 @@ class DDScalar { friend Type operator*(const Scalar a, const Type &b) { return b * a; } - Type &operator*=(const Type &b) { + HYPERJET_INLINE Type &operator*=(const Type &b) { check_equal_size(size(), b.size()); // aliased operands would be read again after being overwritten below @@ -994,7 +1068,7 @@ class DDScalar { return *this; } - Type &operator*=(const Scalar &b) { + HYPERJET_INLINE Type &operator*=(const Scalar &b) { for (index i = 0; i < length(m_data); i++) { m_data[i] = m_data[i] * b; } @@ -1004,7 +1078,7 @@ class DDScalar { // --- div - Type operator/(const Type &b) const { + HYPERJET_INLINE Type operator/(const Type &b) const { using std::pow; check_equal_size(size(), b.size()); @@ -1020,12 +1094,14 @@ class DDScalar { Type result = Type::empty(size()); - binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); + binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } - Type operator/(const Scalar b) const { return Scalar(1) / b * (*this); } + HYPERJET_INLINE Type operator/(const Scalar b) const { + return Scalar(1) / b * (*this); + } friend Type operator/(const Scalar a, const Type &b) { using std::pow; @@ -1036,12 +1112,12 @@ class DDScalar { const auto db = -a / pow(b.f(), 2); const auto dbb = 2 * a / pow(b.f(), 3); - b.unary(b.m_data, f, db, dbb, result.m_data); + b.unary(b.m_data, f, db, dbb, result.m_data); return result; } - Type &operator/=(const Type &b) { + HYPERJET_INLINE Type &operator/=(const Type &b) { using std::pow; check_equal_size(size(), b.size()); @@ -1060,12 +1136,12 @@ class DDScalar { const auto dab = -1 / pow(b.f(), 2); const auto dbb = 2 * a_m_data[0] / pow(b.f(), 3); - binary(a_m_data, b.m_data, f, da, db, daa, dab, dbb, m_data); + binary(a_m_data, b.m_data, f, da, db, daa, dab, dbb, m_data); return *this; } - Type &operator/=(const Scalar &b) { + HYPERJET_INLINE Type &operator/=(const Scalar &b) { operator*=(1 / b); return *this; @@ -1073,7 +1149,7 @@ class DDScalar { // --- arithmetic operations - Type pow(const double b) const { + HYPERJET_INLINE Type pow(const double b) const { using std::pow; Type result = Type::empty(size()); @@ -1082,12 +1158,12 @@ class DDScalar { const auto da = b * pow(this->f(), b - 1); const auto daa = (b - 1) * b * pow(this->f(), b - 2); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type sqrt() const { + HYPERJET_INLINE Type sqrt() const { using std::pow; using std::sqrt; @@ -1097,12 +1173,12 @@ class DDScalar { const auto da = 1 / (2 * f); const auto daa = -da / (2 * this->f()); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type cbrt() const { + HYPERJET_INLINE Type cbrt() const { using std::cbrt; Type result = Type::empty(size()); @@ -1111,26 +1187,26 @@ class DDScalar { const auto da = 1 / (3 * f * f); const auto daa = -da * 2 / (3 * this->f()); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type reciprocal() const { + HYPERJET_INLINE Type reciprocal() const { Type result = Type::empty(size()); const auto f = 1 / this->f(); const auto da = -f * f; const auto daa = -2 * f * da; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } // --- trigonometric functions - Type cos() const { + HYPERJET_INLINE Type cos() const { using std::cos; using std::sin; @@ -1140,12 +1216,12 @@ class DDScalar { const auto da = -sin(this->f()); const auto daa = -cos(this->f()); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type sin() const { + HYPERJET_INLINE Type sin() const { using std::cos; using std::sin; @@ -1155,12 +1231,12 @@ class DDScalar { const auto da = cos(this->f()); const auto daa = -sin(this->f()); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type tan() const { + HYPERJET_INLINE Type tan() const { using std::tan; Type result = Type::empty(size()); @@ -1169,12 +1245,12 @@ class DDScalar { const auto da = f * f + 1; const auto daa = da * 2 * f; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type acos() const { + HYPERJET_INLINE Type acos() const { using std::acos; using std::sqrt; @@ -1186,12 +1262,12 @@ class DDScalar { const auto da = -1 / sqrt(tmp); const auto daa = da * this->f() / tmp; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type asin() const { + HYPERJET_INLINE Type asin() const { using std::asin; using std::sqrt; @@ -1203,12 +1279,12 @@ class DDScalar { const auto da = 1 / sqrt(tmp); const auto daa = da * this->f() / tmp; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type atan() const { + HYPERJET_INLINE Type atan() const { using std::atan; Type result = Type::empty(size()); @@ -1217,12 +1293,12 @@ class DDScalar { const auto da = 1 / (this->f() * this->f() + 1); const auto daa = -da * da * 2 * this->f(); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type atan2(const Type &b) const { + HYPERJET_INLINE Type atan2(const Type &b) const { using std::atan2; Type result = Type::empty(size()); @@ -1236,7 +1312,7 @@ class DDScalar { const auto dab = db * db - da * da; const auto dbb = -daa; - binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); + binary(m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } @@ -1256,8 +1332,7 @@ class DDScalar { const auto dab = -a.f() * b.f() / f3; const auto dbb = a.f() * a.f() / f3; - a.binary(a.m_data, b.m_data, f, da, db, daa, dab, dbb, - result.m_data); + a.binary(a.m_data, b.m_data, f, da, db, daa, dab, dbb, result.m_data); return result; } @@ -1284,15 +1359,15 @@ class DDScalar { const auto dbc = -(b.f() * c.f()) / f3; const auto dcc = (a2 + b2) / f3; - a.ternary(a.m_data, b.m_data, c.m_data, f, da, db, dc, daa, dab, dac, - dbb, dbc, dcc, result.m_data); + a.ternary(a.m_data, b.m_data, c.m_data, f, da, db, dc, daa, dab, dac, dbb, + dbc, dcc, result.m_data); return result; } // --- hyperbolic functions - Type cosh() const { + HYPERJET_INLINE Type cosh() const { using std::cosh; using std::sinh; @@ -1302,12 +1377,12 @@ class DDScalar { const auto da = sinh(this->f()); const auto daa = f; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type sinh() const { + HYPERJET_INLINE Type sinh() const { using std::cosh; using std::sinh; @@ -1317,12 +1392,12 @@ class DDScalar { const auto da = cosh(this->f()); const auto daa = f; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type tanh() const { + HYPERJET_INLINE Type tanh() const { using std::tanh; Type result = Type::empty(size()); @@ -1331,12 +1406,12 @@ class DDScalar { const auto da = 1 - f * f; const auto daa = -2 * f * da; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type acosh() const { + HYPERJET_INLINE Type acosh() const { using std::acosh; using std::sqrt; @@ -1346,12 +1421,12 @@ class DDScalar { const auto da = 1 / (sqrt(this->f() - 1) * sqrt(this->f() + 1)); const auto daa = -da * this->f() / ((this->f() - 1) * (this->f() + 1)); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type asinh() const { + HYPERJET_INLINE Type asinh() const { using std::asinh; using std::sqrt; @@ -1361,12 +1436,12 @@ class DDScalar { const auto da = 1 / sqrt(1 + this->f() * this->f()); const auto daa = -da * this->f() / (1 + this->f() * this->f()); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type atanh() const { + HYPERJET_INLINE Type atanh() const { using std::atanh; using std::pow; @@ -1376,14 +1451,14 @@ class DDScalar { const auto da = 1 / (1 - this->f() * this->f()); const auto daa = 2 * this->f() / pow(this->f() * this->f() - 1, 2); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } // exponents and logarithms - Type exp() const { + HYPERJET_INLINE Type exp() const { using std::exp; Type result = Type::empty(size()); @@ -1392,12 +1467,12 @@ class DDScalar { const auto da = f; const auto daa = f; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type log() const { + HYPERJET_INLINE Type log() const { using std::log; Type result = Type::empty(size()); @@ -1406,12 +1481,12 @@ class DDScalar { const auto da = 1 / this->f(); const auto daa = -da * da; - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type log(const TScalar base) const { + HYPERJET_INLINE Type log(const TScalar base) const { using std::log; Type result = Type::empty(size()); @@ -1420,12 +1495,12 @@ class DDScalar { const auto da = 1 / (this->f() * log(base)); const auto daa = -da / this->f(); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type log2() const { + HYPERJET_INLINE Type log2() const { using std::log; using std::log2; @@ -1435,12 +1510,12 @@ class DDScalar { const auto da = 1 / (this->f() * log(2)); const auto daa = -da / this->f(); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } - Type log10() const { + HYPERJET_INLINE Type log10() const { using std::log; using std::log10; @@ -1450,14 +1525,14 @@ class DDScalar { const auto da = 1 / (this->f() * log(10)); const auto daa = -da / this->f(); - unary(m_data, f, da, daa, result.m_data); + unary(m_data, f, da, daa, result.m_data); return result; } // abs - Type abs() const { return f() < 0 ? -(*this) : *this; } + HYPERJET_INLINE Type abs() const { return f() < 0 ? -(*this) : *this; } // comparison diff --git a/justfile b/justfile index a30acae..f95a55a 100644 --- a/justfile +++ b/justfile @@ -1,6 +1,6 @@ # Format all C++ sources in place format: - uv run --only-group dev clang-format -i include/hyperjet/*.h python/src/*.h python/src/*.cpp test/src/*.cpp benchmark/src/*.cpp + uv run --only-group dev clang-format -i $(git ls-files '*.h' '*.hpp' '*.cpp') # Generate the compilation database that clang-tidy needs compile-db: diff --git a/test/src/test_variants.cpp b/test/src/test_variants.cpp index ad36cdd..87a48da 100644 --- a/test/src/test_variants.cpp +++ b/test/src/test_variants.cpp @@ -252,6 +252,33 @@ TEST_CASE_TEMPLATE("variants: factories", T, D1, X1, D2, X2) { CHECK(vars[1].g(1) == doctest::Approx(1.0)); } +TEST_CASE_TEMPLATE("variants: argument-less factories", T, D1, X1, D2, X2) { + // The dynamic variants start out as scalars carrying no variables, while a + // static one is always its own size. Nothing about an arithmetic result + // depends on this, so it needs a case of its own -- the size returned by the + // argument-less factories was wrong once without a single C++ assertion + // noticing. + const index expected = T::is_dynamic() ? 0 : Size; + + CHECK(T::empty().size() == expected); + CHECK(T::zero().size() == expected); + CHECK(T::constant(S).size() == expected); + + const auto z = T::zero(); + + for (index i = 0; i < z.data_length(); i++) { + CHECK(z.data()[i] == doctest::Approx(0.0)); + } + + const auto c = T::constant(S); + + CHECK(c.f() == doctest::Approx(S)); + + for (index i = 1; i < c.data_length(); i++) { + CHECK(c.data()[i] == doctest::Approx(0.0)); + } +} + TEST_CASE_TEMPLATE("variants: neg", T, D1, X1, D2, X2) { check(-make(A), NegA); }