A C++23 library for differentiating expressions. Build a function over named symbols while the program runs — terms looped over a data file, a model read from configuration, an expression typed by a user, a loss that switches shape on a comparison — and ask it for values, gradients, Jacobians and Hessians. One point at a time, or a batch of thousands in a single call, interpreted or compiled to machine code through LLVM. There are Python bindings over the same runtime, and a header-only compile-time API for expressions whose shape is known when you compile.
#include "ddx.hpp"
#include "rt/equation.hpp"
using namespace ddx;
struct Sample { double t, y; };
const std::vector<Sample> data{{0.0, 2.00}, {1.0, 1.21}, {2.0, 0.74}, {3.0, 0.45}};
// Huber loss of the fit m(t) = a·exp(−b·t): quadratic while a residual is
// inside δ, linear once it is past. `select` is a value, not a branch, so the
// derivative follows whichever arm each residual is on.
const double delta = 0.05;
const auto eq = rt::equation([&] {
const auto a = rt::var("a");
const auto b = rt::var("b");
rt::RTExpression<double> loss = 0.0;
for (const auto &s : data) {
const auto r = s.y - a * exp(-b * s.t);
loss += select(abs(r) < delta, r * r, delta * (2.0 * abs(r) - delta));
}
return loss;
});
*eq.evaluate(2.0, 0.5); // 4.1343959887e-05 — every residual inside δ: the RSS
*eq.jacobian(2.0, 0.5); // {-0.0010757424682, 0.0150678475595}
*eq.hessian(2.0, 0.5); // 2 × 2, row-major
*eq.evaluate(1.8, 0.45); // 2.1570396028e-02 — two residuals past δ, on the linear armContents — Requirements · Using it · Building · Equations · Expressions · Points · Values and derivatives · Batches · Errors · Ownership and threads · Compiling · Saving and loading · Reference · Python · Compile-time expressions · Printing
A C++23 compiler and standard library:
| Compiler | Needs |
|---|---|
| GCC | 14 or newer |
| Clang | 20 or newer, over libstdc++ 14+ — libc++ is not supported |
| MSVC | Visual Studio 2022, /std:c++latest |
Clang compiles the accessor spelling ddx uses from 18 but does not advertise it
until 20, and libc++ has no views::enumerate — hence the two floors.
- CMake 3.26+.
- Boost is downloaded and unpacked at configure time; nothing needs to be
installed, and no compiled Boost library is ever linked. Point
DDX_BOOST_INCLUDEDIRat your own headers to use those instead. - GoogleTest and Google Benchmark are fetched at configure time, so a first configure that builds the tests or the benchmarks wants a network.
-DDDX_BUILD_JIT=ONadditionally needs an LLVM 20 installation, pointed at withLLVM_DIR— to build. The library it produces carries LLVM, and loads on a machine that has none.-DDDX_BUILD_PYTHON=ONadditionally needs Python 3.11+ and pybind11; the module imports NumPy 1.23+ and pydantic 2.7+.
In a project that vendors ddx:
add_subdirectory(ddx)
target_link_libraries(my_app PRIVATE ddx::rt)Or from a package manager — both channels carry the interpreted library, without the JIT:
- vcpkg: an overlay port lives in this repository —
vcpkg install ddx --overlay-ports=<ddx checkout>/contrib/vcpkg/ports— then thefind_packagebelow. The port pins the latest release. - NuGet, for MSVC projects: the
ddxpackage carries the headers and x64 Release and Debug binaries; referencing it wires up the include path, the import library and the DLL copy. The binaries require AVX2.
Or against an installed or built copy:
find_package(ddx CONFIG REQUIRED COMPONENTS rt jit)
target_link_libraries(my_app PRIVATE ddx::jit)cmake --install build/release_jit --prefix /opt/ddxA build tree is consumable without installing:
find_package(ddx CONFIG REQUIRED PATHS /path/to/ddx/build/release NO_DEFAULT_PATH)Asking for a component the build does not have fails at configure time and says which option turns it on.
| Target | Is |
|---|---|
ddx::ddx |
the header-only library — expressions, Equation, forward mode |
ddx::rt |
that plus runtime expressions |
ddx::jit |
ddx::rt with the LLVM backend compiled in |
Components for find_package: ddx, util, ops, md, symbolic, dual,
rt, jit.
cmake -S . -B build
cmake --build build
ctest --test-dir build --output-on-failureOr through the presets:
| Preset | Build type | JIT |
|---|---|---|
debug |
Debug | — |
release |
Release | — |
debug_jit |
Debug | yes |
release_jit |
Release | yes |
debug_jit_tsan |
Debug | yes — under ThreadSanitizer |
python |
Release | yes — and the extension module |
python_no_jit |
Release | — the extension module alone |
cmake --preset release_jit
cmake --build --preset release_jit
ctest --preset release_jitThe version is stated once, in CMakeLists.txt's project() line. CMake
writes util/version.hpp from it — DDX_VERSION_MAJOR/MINOR/PATCH,
DDX_VERSION (100200 for 1.2.0), DDX_VERSION_STRING, and the same as
ddx::version_major, ddx::version_number, ddx::version — the Python module
carries it as ddx.__version__, and a release tag vX.Y.Z must name it.
Every Monday that main has moved since the last tag, CI bumps the patch
number, tags, and publishes the wheels; minor and major bumps are made by hand.
The JIT presets point LLVM_DIR at the Debian/Ubuntu llvm-20 layout;
override it on the command line, which wins over the preset:
cmake --preset release_jit -DLLVM_DIR=/opt/llvm-20/lib/cmake/llvm| Option | Default | Meaning |
|---|---|---|
DDX_BUILD_JIT |
OFF |
compile the LLVM backend into the library |
DDX_BUILD_PYTHON |
OFF |
build the pybind11 extension module |
DDX_BUILD_TESTS |
ON |
build the GoogleTest tests |
DDX_BUILD_BENCHMARKS |
ON |
build the benchmark targets |
DDX_SANITIZE |
off |
thread, address or undefined — instrument the build |
DDX_INSTALL |
on if top-level | generate the install and find_package rules |
ENABLE_NATIVE_ARCH |
ON |
-march=native, else x86-64-v3; x86 only, arm64 takes the compiler's default |
DDX_FP_FLAGS |
ON |
-ffp-contract=fast -fno-math-errno |
DDX_BOOST_INCLUDEDIR |
(empty) | use Boost headers from here instead of fetching |
-ffast-math is not used and not recommended: it changes derivative values.
ddx's own targets compile -fno-exceptions. Nothing in the library throws —
failures come back as values — so what you compile your own code with is your
business.
#include "ddx.hpp"
#include "rt/equation.hpp"
using namespace ddx;Link ddx::rt, or ddx::jit for the compiled batch calls.
rt::equation takes a callback, runs it, and hands back the equation. Symbols
are named inside the callback with rt::var(name) — a string, not a template
argument — and the callback returns the expression it built:
const auto eq = rt::equation([] {
const auto x = rt::var("x");
const auto y = rt::var("y");
return exp(x) * sin(y);
});Return a std::array of expressions for a system, f: ℝⁿ → ℝᵐ:
const auto sys = rt::equation([] {
const auto x = rt::var("x");
const auto y = rt::var("y");
return std::array{x * y + sin(x), exp(x) - y * y};
});
*sys.jacobian(1.3, 0.7); // row-major, m × nAnother scalar type is a template argument on both: rt::var<T>(name) inside
rt::equation<T>(...).
rt::var takes the name as a std::string_view, so it need not be a literal —
the symbols can be whatever the program has just read:
std::vector<std::string> names = read_column_headers(file);
const auto eq = rt::equation([&] {
rt::RTExpression<double> acc = 0.0;
for (const std::string &n : names) {
acc += rt::var(n) * rt::var(n);
}
return acc;
});
*eq.symbols(); // whatever the file named, alphabeticallyNaming the same string twice inside one callback gives the same symbol, so a name that recurs across the data is one variable and not several.
Symbols are ordered canonically — alphabetically by name — which is the
order positional points are read in and the order every result is laid out in.
symbols() lists them in that order, arity() counts them:
*eq.arity(); // 2
*eq.symbols(); // {"x", "y"}An equation can be a string rather than a callback. Include
rt/text/equation.hpp for it. One string per function, so a system is a string
per output:
#include "rt/text/equation.hpp"
const auto eq = *rt::equation("sin(x)*y + 3");
const auto sys = *rt::equation("x*x + y*y - 4", "x*y - 1");These answer result<Equation> rather than an Equation, since a string can
fail to parse where a callback cannot — the * above is that check skipped.
Everything after that is the same — the same points, the same derivatives, the
same batch calls. Free identifiers become the symbols, still ordered
alphabetically, so eq.symbols() here is {"x", "y"}, and a
varmap is the natural point for one.
The grammar is Python's arithmetic: + - * /, ** for exponentiation (^ is
not an operator), parentheses, unary signs, decimals and exponent notation. The
callable functions are the ones in the expression table, spelled
the same. -x**2 is -(x**2) and 2**-1 is legal, as in Python.
Comparisons are infix, as in C++ — < <= > >= == != — and bind
looser than arithmetic, so select(x*x < y + 1, x, y) needs no parentheses.
One comparison per expression: a < b < c is refused rather than read, since Python chains it
into a conjunction and C folds it into (a < b) < c, and picking either
silently would be picking a language.
A string that does not parse comes back as an error — bad_syntax,
unknown_function or wrong_argument_count — rather than as an equation:
const auto bad = rt::equation("x ^ 2"); // `^` is not exponentiation
if (!bad) {
std::println("{}", bad.error()); // errc::bad_syntax
}A string that parses but names no symbol is a different thing: it is an
Equation, and that equation is poisoned with no_graph.
rt::RTExpression<T> is the expression type. A bare number converts to one, so
an accumulator starts at 0.0 and a scalar mixes into arithmetic without
wrapping:
rt::RTExpression<double> acc = 0.0;
for (const auto &s : data) {
acc += (s.y - a * exp(-b * s.t)) * 2.0;
}| Kind | Available |
|---|---|
| Arithmetic | + - * / unary - += -= *= /= |
| Unary | sin cos tan exp log log10 sqrt cbrt abs sign asin acos atan sinh cosh tanh asinh acosh atanh erf |
| Binary | pow atan2 hypot max min |
| Comparison | < <= > >= == != — each answers 1.0 or 0.0 |
| Conditional | select(cond, if_true, if_false) |
They are found by argument-dependent lookup, so they work whether or not rt is
in scope. The compound forms are members, and rebind the handle rather than
mutating what other expressions already refer to.
Where a rule has a choice to make, it is made the same way everywhere: abs
differentiates to 0 at zero and sign to 0 throughout, max and min give
each side half at a tie, and pow(a, b) at a == 0 answers the zeros its
constant arms say — a⁰ is 1 and 0ᵇ is 0, so neither partial is the 0·∞
the textbook rule spells.
select is the only conditional, and it is not a branch: both arms are
evaluated and the condition picks one, so a batch keeps every point on the
same instruction path and the kernel emits a blend. The loss at the top of this
page is the usual shape: a residual tested against a threshold, with a different
expression on each side.
const auto capped = rt::equation([&] {
const auto x = rt::var("x");
return select(x < 1.0, x * x, 2.0 * x - 1.0); // C¹ at x = 1
});The derivative is the taken arm's — d select(c, t, f) = select(c, dt, df) —
and the condition is never differentiated, so the symbols it tests get no
partial through that node. At the switch the function is whatever the two arms
make it: select does not smooth anything, and a discontinuous pair gives a
discontinuous derivative.
A condition is any nonzero value, as in C, and comparisons are ordinary
expressions rather than a separate boolean type — which is why they answer
1.0 and 0.0 and why select(x, …) is legal. Two consequences follow from
IEEE and are worth stating: a comparison against NaN is false (so
NaN == NaN is 0.0), while a NaN used directly as a condition is
true and takes the first arm, exactly as if (nan) does in C.
A comparison answers an expression and never a bool, so if (x < 1.0) does
not compile — but a < b < c does, folded to (a < b) < c as C reads it. A
range test is (a < b) * (b < c).
Only < and <= are operations; the other four are those two read the other
way round, so the graph carries two comparison opcodes and not six.
Every call that needs numbers takes them in five interchangeable spellings:
*eq.jacobian(1.3, 0.7); // positional
*eq.jacobian(std::array{1.3, 0.7}); // any range
*eq.jacobian(named<"y">(0.7), named<"x">(1.3)); // by name, any order
*eq.jacobian("y"_s = 0.7, "x"_s = 1.3); // by name, assignment spelling
*eq.jacobian(std::map<std::string, double>{{"y", 0.7}, {"x", 1.3}}); // a varmapPositional order is alphabetical by symbol name, not the order the symbols appear in the expression.
The named<"y"> and "y"_s spellings key off a name you write down, which a
program that read its symbols out of a file does not have. The varmap is the
spelling for that case: any range of (name, value), the name a run-time string,
matching rt::var(name) on the way in.
std::map<std::string, double> at;
for (const auto &[name, column] : table) {
at[name] = column.back();
}
*eq.evaluate(at);
*eq.jacobian(at);A std::map, an unordered_map, a vector of pairs — anything whose elements
have a .first convertible to std::string_view and a numeric .second. Order
does not matter; each value lands in the slot its name names.
A named point must be whole. Every spelling of one — varmap, named<> and
"x"_s alike — refuses a point that leaves a symbol unreached, with
errc::short_point. Naming a point is not the same as giving one, and a symbol
nothing named would otherwise be read as zero, which is a value rather than a
refusal. A name matching no symbol is errc::unknown_symbol.
point(args…) builds the point vector on its own, in canonical order, so a
point assembled once can be reused:
const auto at = *eq.point(named<"y">(0.7), named<"x">(1.3));
*eq.evaluate(at);Every per-point call answers result<T> — std::expected<T, ddx::error>. The
* in these examples is that check skipped.
| Call | Answers |
|---|---|
evaluate(point) |
result<T>, or result<std::vector<T>> for a system |
jacobian(point) |
result<std::vector<T>>, row-major m × n — so n long when m == 1, which is the gradient |
gradient(point) |
the same block from a graph with no value in it: what only the value needs is not computed |
hessian(point) |
result<std::vector<T>>, dense row-major m × n × n. Each row is its own sweep, so H[i*n+j] and H[j*n+i] can differ in the last ULP — symmetrise before a solver that checks |
univariate_derivative<K>(x0) |
result<T> — the K-th derivative, one symbol and one output only |
const auto g = *eq.jacobian(2.0, 0.5); // {∂f/∂a, ∂f/∂b}
const auto H = *eq.hessian(2.0, 0.5); // H[i * n + j]
H[0]; // 3.1060035856
H[1]; // -3.1441109272For arity()
symbols in symbols() order and output_dim functions:
| Call | Computes | Shape |
|---|---|---|
evaluate(x) |
||
jacobian(x) |
|
|
gradient(x) |
|
|
hessian(x) |
|
|
univariate_derivative<K>(x₀) |
|
scalar |
jvp(v, x) |
|
|
vjp(w, x) |
|
|
hvp(v, x) |
|
The last three are the matrix-free products. Each is one number per output, and none of them ever forms the matrix it is named after.
const std::vector<double> v{1.0, 0.0};
*eq.jvp(v, 2.0, 0.5); // J·v — m long: how f moves if x moves along v
*eq.vjp(v, 2.0, 0.5); // wᵀJ — n long: ∇(w·f), one weight per function
*eq.hvp(v, 2.0, 0.5); // H·v — n long: how ∇f moves if x moves along vhvp also takes a symbol name, for the common case where eq.hvp("a", x) is a. It refuses with
errc::unknown_symbol if no such symbol exists.
*eq.hvp("a", 2.0, 0.5); // {∂²f/∂a², ∂²f/∂b∂a}Why they are worth asking for, and when they are not. Each answers in one
number per output and needs no matrix: hvp is hessian() is vjp is jacobian() is one column
per structural nonzero; jvp is
What that buys is storage and bandwidth, not time. Measured
(benchmarks/jit/benchmark_hvp.cpp), one hvp costs about the same as one
whole hessian — within 10% across
-
Want the whole matrix? Use
hessian(). Asking for$n$ products to assemble it is far slower than one colouring — at$n = 64$ , roughly 70×. -
Want one product, or cannot afford
$n^2$ storage? Usehvp. At$n = 64$ the Hessian block is 4096 columns per point against 64.
A direction is one value per symbol (errc::wrong_direction and computes nothing.
const auto in_b = rt::equation([] { return sin(rt::var("b")); });
*in_b.univariate_derivative<4>(0.5); // sin⁗(0.5)univariate_derivative sweeps the graph once in truncated-Taylor arithmetic,
and every operation above is defined there — max, min, abs, select and
the comparisons included, taking the side the point is on. Should a scalar type
ever lack one of them, the call answers errc::unsupported_scalar rather than a
zero.
Values, gradient and Hessian are prepared separately, each the first time it is asked for. A caller who only ever evaluates never pays for a gradient.
Hand the factory a cache and an equation keeps its last call per kind. Ask again at the same point and the answer comes back off it; move one symbol and only the part of the graph that symbol reaches is computed again. Nothing else changes: the numbers are the ones the sweep gives, bit for bit, and there is no call to make -- a point already answered is simply answered.
const auto eq = ddx::rt::equation(model); // off
const auto eq = ddx::rt::equation(model, ddx::rt::LastValue{}); // onIt is off by default, and it earns its place where points repeat or arrive one
coordinate at a time -- a line search, a finite difference, a minimiser asking
for the value and then the gradient at the same x:
const auto f = *eq.evaluate(x); // swept
const auto g = *eq.gradient(x); // its own lane, so its own first call
const auto h = *eq.evaluate(x); // nothing swept at allEvery spelling that takes a point is covered; a batch is not, being amortised already. A lane a kernel answers keeps the point but not the tape, so it serves a repeat and computes a moved point whole -- compiled code has no way to run part of itself.
LastValue holds one slot per kind behind a shared lock: any number of threads
can be served the same point at once, and one that arrives while another is
writing sweeps for itself rather than waiting. The object is copied in, so it
may be a temporary.
To supply your own, model ddx::rt::CValueCache: answer active(), and hand out
a read lease and a write lease per Want, each lending the point, the output
values and the tape of one remembered call. active() is asked before anything
else and is how a cache says no call of yours can hit -- a false there costs
the equation nothing at all, not even reading the point. ddx::rt::Extent says
what a slot has to hold and which frozen graph it holds it for; two calls belong
to the same entry only where their whole Extent matches, so a re-freeze or a
grown arena parts them without the cache having to know why.
The batch calls take columns: one input column per symbol, one output pointer
per output column, every column n points long. They fill the outputs in place
and answer result<void>.
result<void> evaluate(xs, values, n);
result<void> jacobian(xs, values, partials, n);
result<void> hessian (xs, values, partials, hessians, n);
// The seeded products take a second block of input columns: `vs` is the
// direction at each point, one column per symbol, and `ws` one per function.
result<void> jvp(xs, vs, values, products, n); // J v, m columns
result<void> vjp(xs, ws, values, products, n); // wᵀJ, n columns
result<void> hvp(xs, vs, values, partials, products, n); // H v, n columnshvp fills partials with
Ask the equation how many output columns each block wants, and size the buffers by that:
const std::vector<double> as{1.8, 1.9, 2.0, 2.1}; // candidate a values
const std::vector<double> bs{0.45, 0.48, 0.50, 0.52}; // candidate b values
const std::size_t n = as.size();
std::vector<double> f(n * *eq.value_columns());
std::vector<double> j(n * *eq.jacobian_columns());
const std::vector<const double *> xs{as.data(), bs.data()};
const std::vector<double *> values{f.data()};
const std::vector<double *> partials{j.data(), j.data() + n};
if (const auto ok = eq.jacobian(xs, values, partials, n); !ok) {
std::println("{}", ok.error());
}
f[2]; // 4.1343959887e-05 — the loss at (2.0, 0.50)
j[2]; // ∂loss/∂a there
j[n + 2]; // ∂loss/∂b there
f[0]; // 2.1570396028e-02 — (1.8, 0.45): two terms on the linear arm
j[0]; // -0.1617863027 — and the gradient is theirs, ±2δ·∂r/∂aEvery point of a batch runs the same instructions: the select in the model
computes both arms and keeps one per point, so a batch that straddles δ costs no
more than one that does not, and compiles to a single blend.
A column count that does not match answers errc::wrong_column_count, and
nothing is written.
Both derivative blocks are sparse, so jacobian_columns() and
hessian_columns() are what size them — asking costs nothing — and a pattern
says which cell each column is. A caller who wants the dense matrix uses the
per-point calls above, which undo both compressions on the way out.
The Jacobian keeps only the cells that structurally exist. ∂fᵢ/∂xⱼ gets a
column when the derivative is something other than the literal zero — a symbol
a function does not mention has no column, and neither has one whose partial
folded away. jacobian_pattern() places them:
const rt::Sparsity &pattern = eq.jacobian_pattern()->get();
pattern.nonzeros(); // == *eq.jacobian_columns()
for (const rt::Cell &c : pattern.entries()) { // every cell that exists
use(c.row, c.column, j[c.slot * n + point]); // ∂f(row)/∂x(column)
}
const double dfi_dxj = // at() is nullopt off the pattern
pattern.at(i, j)
.transform([&](std::size_t cell) { return j[cell * n + point]; })
.value_or(0.0);For a system whose functions each touch a few symbols that is most of the matrix: a tridiagonal residual over six variables has 16 columns rather than 36, and the graph, the kernel and the compile shrink with it.
The Hessian takes a fourth block, compressed by colour rather than by pattern —
one output only. hessian_cell(i, j) places those:
std::vector<double> h(n * *eq.hessian_columns());
const auto hessians =
std::views::iota(0uz, *eq.hessian_columns()) |
std::views::transform([&](std::size_t k) { return h.data() + k * n; }) |
std::ranges::to<std::vector>();
const auto ok = eq.hessian(xs, values, partials, hessians, n);
const auto cell = eq.hessian_cell(i, j); // std::optional
const double d2 = cell ? h[*cell * n + point] : 0.0;evaluate is not the Jacobian call with the partials thrown away: columns
nobody asks for are work nobody does.
rt::equation always hands back an Equation. When the callback or the string
could not produce one, that equation is poisoned: it carries the error rather
than a function, and every call on it answers with that error.
const auto eq = rt::equation([] { return rt::var("x") * 2.0; });
if (const auto bad = eq.status()) {
std::println("{}", *bad); // errc::no_arena
return;
}status() is a std::optional<ddx::error>; poisoned() is the same question
without the reason.
Calls answer result<T> on top of that: the symbol list exists only at run
time, so an arity mismatch or an unknown name is a genuine runtime failure.
const auto j = eq.jacobian(1.0, 2.0); // two values, one symbol
if (!j) {
std::println("{}", j.error()); // errc::wrong_arity
}errc |
Means |
|---|---|
wrong_arity |
the point does not supply one value per symbol |
wrong_direction |
a direction does not supply one value per symbol, or a covector one per function |
short_point |
a range point is shorter than the symbol list, or a named point leaves a symbol unreached |
unknown_symbol |
a named point uses a name the equation does not have |
index_out_of_range |
the index does not name a symbol of this equation |
wrong_column_count |
a batch block has the wrong number of columns |
no_arena |
a symbol was named outside an rt::equation callback |
no_graph |
the expression is a bare literal, naming no function |
sealed_arena |
the symbols already back an equation, so they are final |
not_univariate |
univariate_derivative on more than one symbol |
unsupported_scalar |
univariate_derivative over an operation its Taylor arithmetic does not define |
bad_syntax |
a text equation the grammar does not accept |
unknown_function |
a text equation calls a function that does not exist |
wrong_argument_count |
a text equation calls one with the wrong arity |
archive_io |
the file could not be read or written |
bad_archive |
not a ddx file, or a format this build does not read |
archive_corrupt |
the file's checksum or structure does not hold |
archive_mismatch |
the file loads, but does not describe this equation |
jit_target, jit_module, jit_object, jit_verify, jit_lookup |
the compiler could not produce or link a kernel |
An errc and an error both format and stream as their text, so
std::println("{}", j.error()) is the whole of reporting one. Nothing here
throws.
The accessors that answer no result answer std::optional instead —
arity(), symbols(), value_columns(), jacobian_columns(),
hessian_columns(), hessian_colors() and jacobian_pattern() are nullopt
on a poisoned equation.
A bare count could not say it: an equation over a literal-only expression
legitimately has no symbols and no output columns.
The equation owns everything it needs. It is move-only: return it from a function, store it in a class, keep it as long as you like.
auto make_model(const std::vector<Sample> &data) {
return rt::equation([&] {
const auto a = rt::var("a");
const auto b = rt::var("b");
rt::RTExpression<double> rss = 0.0;
for (const auto &s : data) {
const auto r = s.y - a * exp(-b * s.t);
rss += r * r;
}
return rss;
});
}
const auto eq = make_model(load()); // `data` is gone; `eq` is completeSymbols belong to the callback. rt::var is meaningful while the callback
that named it runs. Build the expression there and return it; do not store an
RTExpression and use it afterwards. A symbol named outside a callback gives a
poisoned equation carrying errc::no_arena.
Batch columns are yours. The equation writes through the pointers you give
it and keeps none of them; every buffer must be alive for the duration of the
call, and each column must hold at least n elements.
Building is independent per thread: two threads may run rt::equation callbacks
at the same time. An equation prepares itself on the first call of each kind —
make that first call before handing one to several threads, or give each thread
its own. A batch call splits across threads by slicing the columns, each thread
getting its own offset pointers and its own n:
const std::size_t chunk = (n + threads - 1) / threads;
std::vector<std::jthread> pool;
for (std::size_t t = 0; t < n; t += chunk) {
const std::size_t m = std::min(chunk, n - t);
pool.emplace_back([&, t, m] {
const std::vector<const double *> xs{as.data() + t, bs.data() + t};
const std::vector<double *> values{f.data() + t};
const std::vector<double *> partials{j.data() + t, j.data() + n + t};
(void)eq.jacobian(xs, values, partials, m);
});
}Configure with -DDDX_BUILD_JIT=ON and link ddx::jit, and the batch calls can
run compiled code. The spellings do not change, the answers do not change, and
there is no compiler or kernel object for you to hold.
Nothing compiles unless you ask. An equation left alone interprets, which costs no compiler and runs within ~1.4x of a kernel. Asking is what starts the build, so it overlaps whatever you do next, and no call ever waits for it:
auto eq = rt::equation([] {
const auto x = rt::var("x");
const auto y = rt::var("y");
return x * log(x) + y * exp(x * y);
});
eq.options({.backend = rt::Backend::Compile, .points = 4096});
// Starts the compile here and now. Returns immediately; nothing is ready yet.
eq.jacobian(xs, f, g, n); // compile in flight -> interpreted
// ... the kernel lands about here ...
eq.jacobian(xs, f, g, n); // compiled, and from here on
eq.uses_kernel(); // trueValues, gradient and Hessian compile separately, each launched the moment it is first needed — so a caller who only evaluates never compiles a gradient.
Results either side of the switchover agree to the bit, with one exception worth knowing. The compiled path sums a reduction spine of sixteen terms or more in blocks, where the interpreted one sums it left to right — k dependent adds are a latency no lane width hides. Addition is associative in exact arithmetic and not in floating point, so a result reached through such a spine can move in its last bits the moment the kernel lands: measured at 2 ULP over forty terms and 4 over eighty. Shorter spines are not rewritten and cannot move, and the Hessian lane is built from one graph either way, so it never moves.
wait_for_kernel() is the only call that blocks, and only because it was asked
to:
eq.wait_for_kernel(); // true once it lands
eq.jacobian(xs, f, g, n); // compiled, guaranteedSetting Options to the value it already holds does nothing. Changing it
discards the compiled code and abandons a compile still in flight rather than
waiting for it — including a change back to Backend::Interpret, after which no
compiler is asked at all.
Backend |
Is |
|---|---|
Interpret |
the default. No compiler is asked, so a program that never says otherwise never loads LLVM |
Compile |
start compiling now |
Adapt |
start compiling once the batch traffic has paid for it — for a caller who cannot say up front which of their equations are hot |
Either compiling backend answers first with a quick kernel and replaces it with
a better one; every level agrees to the bit, so the swap is invisible.
kernel_level() says which one is answering, and under Adapt, warming()
answers how far along the next step is:
eq.options({.backend = rt::Backend::Adapt});
if (const auto w = eq.warming()) {
std::println("{} of {} points", w->points, w->threshold);
}
*eq.kernel_level(); // 0 at first, then the level you asked for| Field | Default | Is |
|---|---|---|
backend |
Interpret |
Interpret, Compile or Adapt |
points |
1 |
the batch you intend to hand one call — stated, since the kernel is built before any call exists to infer it from |
codegen.lanes |
Lanes::derived() |
points per loop iteration; derived is the host's register width, scalar for a batch too short to fill one. Lanes::scalar() or *Lanes::exactly(w) states one. Every width gives the same bits |
codegen.opt_level |
follows the build type | LLVM's IR pipeline, Level::O0 to Level::O3 — O3 in a Release build, O1 in a Debug one |
codegen.codegen_level |
Level::O1 |
LLVM's codegen, Level::O0 to Level::O3 — the knob that trades kernel speed for compile time |
codegen.slp |
false |
pack independent subexpressions within one point; model-dependent, hence off |
codegen.loop_vectorize |
false |
loop vectorisation, on a loop already emitted lanes wide |
codegen.veclib |
None |
vector math library for transcendentals; Libmvec trades ~0.5 ULP for ~4, and a derived lanes is then the widest it serves (four) |
codegen.contract |
follows DDX_FP_FLAGS |
fold a multiply feeding an add into one rounding |
warm_points |
65536 |
batch points that buy the first step, under Adapt |
hot_points |
1048576 |
further points that buy the top one, under Adapt |
retain_object |
true |
keep the compiled object so save can write it |
cache_dir |
(empty) | keep compiled objects here between runs; a second run links instead of compiling |
time_passes |
false |
per-pass timing to stderr |
codegen is everything the emitter reads, and so the identity a stored kernel
is matched against and the object cache is keyed on; the fields around it are
policy and change no machine code.
points decides the lane width and nothing else: a call carrying some other
number is answered correctly, just not by the kernel that number would have
built. At one point a scalar kernel is 1.2x to 4.1x quicker than the host's
vector width, and over a batch the reverse by 3x.
eq.options({.backend = rt::Backend::Compile, .points = 1}); // a gradient per step
eq.options({.backend = rt::Backend::Compile, .points = 4096}); // a batch at a timeNaming a cache_dir is the cheapest thing you can do about compile time: a run
that finds its object there links it instead of compiling, which is roughly
three orders of magnitude quicker. Failing that, codegen.codegen_level — not
codegen.opt_level — is the knob that moves a compile time.
eq.options({.backend = rt::Backend::Compile, .points = 4096, .cache_dir = "/var/cache/ddx"});Preparing an equation for derivatives is work you can do once. save writes it
to a file, and load reads it back without a model:
const auto eq = rt::equation([] { return exp(rt::var("x")) * rt::var("y"); });
eq.save("f.ddx");
const auto same = *rt::load("f.ddx"); // no model, nothing rebuilt
const auto sys = *rt::load<double, 3>("s.ddx"); // three functionsThe output count is a template parameter because it is part of the type, and a file holding some other number is refused rather than adapted to.
Pair a model with a file and it becomes a cache — the lambda is the model as you would write it anyway, and the path is only where the result is kept:
const auto eq = rt::equation("f.ddx", [] {
const auto x = rt::var("x");
return exp(x) * x;
});The first run builds and writes the file. A later one finds the file still
describes the model and takes the work off disk. Editing the model rebuilds and
overwrites rather than trusting a stale file, and it never refuses: an absent
file is the first run, a stale one is a rebuild. loaded() says which happened.
Two questions about a file, and they are different:
rt::verify("f.ddx"); // result<void> — can this build read it at all?
eq.verify("f.ddx"); // result<void> — ...and is it *this* equation?Files are binary and little-endian, checked before anything in them is used, so one written by a different version of ddx is refused rather than misread.
Compiled code travels in the file too:
eq.options({.backend = rt::Backend::Compile});
eq.wait_for_kernel();
eq.save("f.ddx"); // the kernel goes with it
const auto warm = *rt::load("f.ddx");
warm.uses_kernel(); // true, with nothing having compiledMeasured at 64, 128 and 256 variables: 38/79/179 ms to build and compile against
0.40/0.63/1.12 ms to load. retain_object is on by default so the code is there
when a save comes; a caller that never saves turns it off, and the file then
carries the equation and no code.
Stored code runs only where the equation, the host and the options it was built under all still agree. Anything else compiles instead — a mismatch is answered by compiling, never by running the wrong code.
| Member | Answers |
|---|---|
save(path) |
result<void> — write this equation |
verify(path) |
result<void> — does path hold this equation? |
Equation::load(path) |
result<Equation> — read one |
loaded() |
whether this equation came off disk |
rt::verify<T>(path) |
result<void> — is path readable by this build? |
Five functions, a gradient with respect to every variable at one point, the
point changing every call. Every cell is the time for the gradient divided by
the time to evaluate the same function once at plain double — the overhead
factor, so lower is better and below 1.00x means the derivatives came cheaper
than the function — as the median of 3 interleaved repetitions, with the
spread (max − min) / median beside it. The last column is that one evaluation.
| arm | what it is |
|---|---|
ddx-rt/r |
Equation::gradient, the graph interpreted |
ddx-jit/r |
the same graph compiled through LLVM, one point wide |
casadi/r |
CasADi's SX gradient Function, interpreted, through Function.buffer() |
casadi/j |
the same Function with jit=True, its generated C compiled by the system compiler |
adept/r |
Adept 2.1.3, new_recording and compute_adjoint per point |
Regular solid solution
| n | ddx-rt/r | ddx-jit/r | casadi/r | casadi/j | adept/r | one evaluation, ns |
|---|---|---|---|---|---|---|
| 16 | 3.61x ±5% | 1.03x ±2% | 11.59x ±7% | 5.79x ±6% | 9.22x ±11% | 195 |
| 32 | 3.94x ±3% | 1.02x ±1% | 10.70x ±11% | 2.48x ±6% | 10.19x ±3% | 650 |
UNIQUAC
| n | ddx-rt/r | ddx-jit/r | casadi/r | casadi/j | adept/r | one evaluation, ns |
|---|---|---|---|---|---|---|
| 16 | 4.57x ±2% | 0.92x ±1% | 8.96x ±2% | 1.92x ±4% | 6.33x ±5% | 807 |
| 32 | 4.43x ±3% | 0.93x ±1% | 7.85x ±3% | 1.55x ±1% | 6.90x ±1% | 2250 |
Peng-Robinson — n counts N−1 mole fractions plus Z.
| n | ddx-rt/r | ddx-jit/r | casadi/r | casadi/j | adept/r | one evaluation, ns |
|---|---|---|---|---|---|---|
| 17 | 3.10x ±5% | 0.48x ±2% | 9.28x ±3% | 1.69x ±3% | 4.42x ±3% | 857 |
| 33 | 2.89x ±6% | 0.45x ±2% | 7.84x ±17% | 1.03x ±3% | 4.18x ±1% | 3340 |
Mixed Solvent Electrolyte
| n | ddx-rt/r | ddx-jit/r | casadi/r | casadi/j | adept/r | one evaluation, ns |
|---|---|---|---|---|---|---|
| 16 | 5.16x ±3% | 0.95x ±1% | 12.94x ±1% | 1.98x ±2% | 7.96x ±2% | 1164 |
| 32 | 5.20x ±3% | 1.03x ±2% | 11.33x ±1% | 1.55x ±2% | 8.69x ±2% | 3589 |
The op-coverage function
| n | ddx-rt/r | ddx-jit/r | casadi/r | casadi/j | adept/r | one evaluation, ns |
|---|---|---|---|---|---|---|
| 16 | 2.62x ±3% | 0.71x ±1% | 2.96x ±3% | 1.39x ±5% | 2.98x ±6% | 1616 |
compare/ is the harness — the functions, the arms, the gate that checks every
arm against every other before anything is timed, and the synthetic families
behind --trend that the interpreter is tuned against.
ddx::rt::equation(callback) → Equation. All of it is const and
thread-safe except options().
| Member | Answers |
|---|---|
poisoned(), status()
|
whether the build failed, and with what |
arity() |
std::optional<std::size_t> — symbol count, n
|
symbols() |
std::optional<std::span<const std::string>> — canonical order |
point(args…) |
result<std::vector<T>> — a point in canonical order, from any of the five spellings |
evaluate(point) |
result<T>, or result<std::vector<T>> for a system |
jacobian(point) |
result<std::vector<T>>, dense row-major m × n |
gradient(point) |
the same, from a graph with no value block |
hessian(point) |
result<std::vector<T>>, dense row-major m × n × n |
univariate_derivative<K>(x0) |
result<T>, one symbol and one output only |
jvp(v, point) |
result<std::vector<T>>, |
vjp(w, point) |
result<std::vector<T>>, |
hvp(v, point) |
result<std::vector<T>>, |
hvp(name, point) |
the same, along the basis vector for that symbol |
evaluate(xs, f, n) |
result<void> — a batch of n
|
jacobian(xs, f, g, n) |
result<void> |
gradient(xs, g, n) |
result<void> |
hessian(xs, f, g, h, n) |
result<void>, one output only |
jvp(xs, vs, f, p, n) |
result<void> |
vjp(xs, ws, f, p, n) |
result<void> |
hvp(xs, vs, f, g, p, n) |
result<void>, one output only — g gets |
value_columns(), jacobian_columns(), hessian_columns()
|
std::optional<std::size_t> — what sizes the batch buffers |
jvp_columns(), vjp_columns(), hvp_columns()
|
the same, for the seeded products — m, n, n
|
jacobian_pattern() |
std::optional<std::reference_wrapper<const rt::Sparsity>> — which (function, symbol) cell each Jacobian column is |
hessian_cell(i, j) |
std::optional<std::size_t> — which Hessian column holds H(i, j), or none |
hessian_colors() |
std::optional<std::size_t> — groups in the Hessian's compression |
options(opts), options()
|
set or read the compile options; setting returns *this
|
uses_kernel(), kernel_level()
|
whether a batch call runs compiled code, and at which level |
warming() |
under Adapt, points seen against the next threshold |
wait_for_kernel() |
block until a compile in flight has landed |
save(path), verify(path)
|
write this equation; ask whether a file holds it |
load(path), loaded()
|
read one; whether this one was read |
The same runtime, as an extension module. A release carries one wheel per platform and CPython version (3.11–3.14) and a source distribution. A wheel needs nothing installed beside it — LLVM is inside the library:
uv pip install https://github.com/reach2sayan/ddx/releases/download/v1.2.1/ddx-1.2.1-cp312-cp312-manylinux_2_28_x86_64.whl| Wheel | JIT |
|---|---|
| Linux x86_64 (glibc 2.28+) | yes |
| Windows x64 | no — calls interpret |
There is no macOS wheel: the tree uses C++23 ranges that libc++ does not have
(views::enumerate, cartesian_product, chunk, stride, ranges::fold),
so on a Mac it builds from source with a libstdc++ toolchain and without the
JIT.
Building from source instead — pip install ., or a preset — needs a C++23
compiler and, for the JIT, LLVM 20's archives with a static zlib and zstd;
scripts/build_llvm.py builds that set from source into a prefix, and is what
the Linux wheel uses:
pip install . # wheel; JIT on, so building needs LLVM 20
cmake --preset python # in-tree, JIT
cmake --preset python_no_jit # in-tree, no LLVMddx.has_jit says whether the copy you have was built with the LLVM backend.
Everything below works either way; without it, calls interpret.
equation takes a model — a callable of no arguments returning one expression,
or a tuple of them for a system — as a call or a decorator:
import ddx
@ddx.equation
def f():
x = ddx.var("x")
y = ddx.var("y")
return ddx.exp(x) * ddx.sin(y)Or a string, or a list of strings for a system, in the same grammar the C++ side reads:
g = ddx.equation("exp(x) * sin(y)")
sys = ddx.equation(["x*x + y*y - 4", "x*y - 1"])Symbols are named inside the model with ddx.var(name), and they order
alphabetically as they do in C++. f.symbols lists them, f.arity counts them,
f.outputs counts the model's outputs. Arithmetic operators work on
Expression, and a bare number mixes in; the free functions are the same set
the C++ side has.
A point is a sequence, or a dict keyed by symbol name. A 2-D array is a
batch — shape (symbols, points), one row per symbol:
f.jacobian([2.0, 3.0]) # positional, alphabetical
f.jacobian({"x": 2.0, "y": 3.0}) # by name
f.jacobian(np.array([[2.0], [3.0]])) # a batch of one
f.jacobian(np.array([[2.0, 2.5], [3.0, 3.5]])) # a batch of twoEach call returns everything up to and including what it names, so a gradient never costs a second evaluation:
| Call | Returns |
|---|---|
evaluate(x) |
f |
jacobian(x) |
(f, J) |
gradient(x) |
J alone, from a graph that does not compute f
|
hessian(x) |
(f, J, H) |
jvp(v, x) |
|
vjp(w, x) |
|
hvp(v, x) |
|
value, gradient = f.jacobian([2.0, 3.0])
value, gradient, hessian = f.hessian([2.0, 3.0])
hessian.shape # (2, 2)
value, gradient, hv = f.hvp([1.0, 0.0], [2.0, 3.0])
hv.shape # (2,) — H·v, without forming HShapes follow the point: a single point gives a scalar f, an (n,) gradient
and an (n, n) Hessian; a batch of p appends that axis, giving (p,),
(n, p) and (n, n, p). A system prepends its output axis. Unlike the C++
batch calls, the Hessian arrives dense — the compression is undone on the
way out, so nothing here needs hessian_columns.
The three products never form the matrix they are named after, which is what
makes them worth having when (symbols, points) batch alongside a batched
point. hvp needs a single-output model, as hessian's compressed path does.
The calls above allocate their answers. buffer(x) binds the point and the
answers once and hands back a Call: write the next point into x, call it,
read the blocks back. Same arrays every time, so nothing is allocated per call.
call = f.buffer(np.array([2.0, 3.0]))
for _ in range(steps):
call() # fills call.value and call.jacobian
call.x[:] = next_point(call.jacobian)want chooses how far it goes, and a block nobody asked for is one nobody
computes:
want |
Fills |
|---|---|
Want.VALUE |
value |
Want.JACOBIAN (default) |
value, jacobian |
Want.GRADIENT |
jacobian |
Want.HESSIAN |
value, jacobian, hessian |
Reading a block the call did not ask for raises errc.wrong_column_count, and
so does binding Want.HESSIAN on a system — a frozen graph carries one
colouring, so m Hessians are read off the arena a point at a time.
Shapes are the ones the allocating calls answer with, value included: one
output at one point is a float, and everything else is an array. The point is
bound as an array whatever was passed, so call.x is writable even when the
argument was a list.
remember=True gives the equation the same last-call cache the C++ side gets
from LastValue, on the same terms: a repeated point is answered off the last
one, a point one symbol away sweeps only what that symbol reaches, and the
numbers are unchanged.
f = ddx.equation(model, remember=True)
value = f(x) # swept
grad = f.gradient(x) # its own lane
again = f(x) # nothing sweptIt applies to a point at a time and not to an array of them, and a lane a kernel answers serves a repeat without amending a moved point.
Python raises where C++ returns result<T>. ddx.Error is a RuntimeError
carrying the same code:
try:
f.jacobian({"z": 1.0})
except ddx.Error as e:
print(e) # "no symbol of that name"
e.code is ddx.errc.unknown_symbol # Trueerrc is an IntEnum, so it compares and formats as its number — e.code.name
is the spelling. The codes are the ones in the table above; the
exception's own message is the text.
Options is a frozen pydantic model of the same fields as jit::Options,
validated on the way in. eq.options reads and assigns it; eq.compile() sets
backend=COMPILE, waits for the kernel, and returns the equation, so a
configure-and-use reads in one line:
f.compile(points=batch.shape[1]).jacobian(batch)
f.uses_kernel # True, once it has landed
f.options = ddx.Options(backend=ddx.Backend.INTERPRET) # discards the kernelcompile() blocks by construction — it is wait_for_kernel() with the options
set first. Assigning options does not: calls interpret until the kernel lands
and switch over when it does, as in C++.
Equations save and load here too, over the same file format — a file written by C++ loads in Python and the other way round:
eq.save("f.ddx")
same = ddx.load("f.ddx") # no model runs, nothing is rebuilt
@ddx.equation # or pair a model with a file, as a cache
def model() -> ddx.Expression:
x, y = ddx.var("x"), ddx.var("y")
return ddx.exp(x) * y
cached = ddx.equation(model, cache="f.ddx")
cached.loaded # False the first run, True afterA text equation caches the same way: ddx.equation("exp(x) * y", cache="f.ddx").
save, load and verify raise ddx.Error rather than answering False:
unreadable, unloadable and "a different equation" are three different errc
values, and only the code says which.
| Member | Is |
|---|---|
arity, outputs, symbols
|
properties — symbol count, output count, canonical names |
evaluate(x), __call__(x)
|
f at the point or batch |
jacobian(x) |
(f, J) |
gradient(x) |
J alone |
jvp(v, x), vjp(w, x), hvp(v, x)
|
|
hessian(x) |
(f, J, H), dense |
options |
property — read or assign an Options
|
compile(**fields) |
set Options, block for the kernel, return self |
uses_kernel, wait_for_kernel(*, want)
|
whether a call runs compiled code, and blocking for it — for the Jacobian lane unless want names another |
hessian_colors |
groups in the Hessian's compression |
buffer(x, *, want) |
a Call bound to its buffers, for a loop |
to_dot(*, all=False) |
the expression in Graphviz form; all=True draws the pruned nodes too |
nodes(*, want) |
how many nodes a call for want evaluates |
save(path), verify(path)
|
write this equation; raise unless path holds it |
loaded |
property — whether this equation was read rather than built |
ddx.load(path) reads one, and ddx.equation(model, cache=path) builds or reads
as the file allows.
The same derivatives are available on expressions whose shape is known when you
compile, header-only and constexpr throughout — #include "ddx.hpp", link
ddx::ddx. Symbols are named by template argument, and Equation carries the
whole API:
constexpr auto x = var<"x">;
constexpr auto y = var<"y">;
constexpr auto eq = Equation{x * y + 2.0 * x};
constexpr auto g = eq.jacobian(std::array{2.0, 3.0});
static_assert(g[0] == 5.0); // ∂f/∂x
static_assert(g[1] == 2.0); // ∂f/∂yPoints come in the same four spellings, and results are std::array or a tensor
rather than result<std::vector<T>> — the shapes are in the type, so nothing
has to be checked at run time.
| Member | Answers |
|---|---|
evaluate(point) |
f at the point |
jacobian(point) |
J, J[i, j] = ∂fᵢ/∂xⱼ; ∇f itself when m == 1 |
hessian(point) |
∇²f — on symbols declared var<"x", dual> |
derivative_tensor<K>(point) |
all K-th order partials, any K ≥ 1 |
univariate_derivative<K>(x0) |
the K-th derivative of a one-symbol function |
constexpr auto H = Equation{var<"a", dual> * var<"b", dual>}.hessian(2.0, 3.0);
H[0, 1]; // 1.0
Equation{x * y}.derivative_tensor<3>(1.0, 2.0); // rank 3record(named<"n">(3), named<"x">(1.5)) collects keyword arguments into a value
with the keys in its type, readable with m.get<"n">() or m["x"_s].
Compile-time expressions and Equations are formattable and streamable:
std::format("{}", x * y + sin(x)); // "x * y + sin(x)"
std::format("{::.3f}", 2.0 * x); // "2.000 * x" — the spec applies to every number
std::cout << (x - y * x) << '\n'; // "x - y * x"An Equation prints each output with its derivative row underneath. A symbol
held constant for one partial prints with a _c suffix:
f0: x * y
grad: y_c, x_c
Slot 0 of an Equation is the function and slot k > 0 is ∂f/∂xₖ in canonical
order, each an expression you can print or store:
auto eq = Equation{x * y};
std::format("{}", eq[idx<0>()]); // "x * y"
std::format("{}", eq[idx<1>()]); // "y_c" — ∂f/∂xBoost Software License 1.0. Suggestions and pull requests are welcome.