odelia architecture: the XAD Tape runtime and how downstream packages link it

Read this before changing how XAD is compiled, how odelia.so/odelia.dll is loaded, or R/zzz.R. These pieces form a contract that packages which LinkingTo: odelia depend on. Breaking it produces link-time or load-time failures in downstream packages (e.g. plant) that do not show up in odelia’s own checks.

For automatic differentiation, read the article Reverse mode: one solve, every parameter (vignettes/articles/reverse-mode.Rmd): what a run records, how the sweep walks it, what a System must provide, and the order to read the headers in. The contract a System implements is the set of concepts in inst/include/odelia/ode_interface.hpp, which the compiler checks. This document is about how the headers fit together and how the Tape runtime compiles and links.

Map of the headers

inst/include/odelia/, in dependency order; a header includes only those above it.

Header	Holds
`ode_util.hpp`	`util::stop` / `util::warning` (R-free), `to_passive`, `DomainError`
`tangent.hpp`	`tangent_scalar<T>`, the forward-mode scalar; refuses a tangent above an adjoint
`adjoint_scalar.hpp`	`adjoint_scalar<T>` and `adjoint_tape<T>`, the reverse-mode scalar and its tape; `visit_active` and `count_active_slots`, the walk over the active values a class holds
`with_slope.hpp`	`with_slope<T>`, a value paired with its slope; no includes
`spline.hpp` → `interpolator.hpp` → `drivers.hpp`	cubic Hermite backend, `hermite_interpolator`, `Drivers` (time-varying forcing)
`ode_control.hpp`	`OdeControl`: tolerances, step bounds, the `gsl` / `hairer` controller
`ode_interface.hpp`	the System concepts (`HasOdeTime`, `Rebindable`, `SolvesForValues`, `ChecksState`, `Sweepable`, `HasSubsystem`, `ScalesTolerances`), `schedule_row`, `step_record`, `be_at_step`, `apply_insertion`
`adjoint.hpp`	`adjoint_batch`, `adjoint_system`, `vector_jacobian_product`, `state_and_parameter_adjoints`: the transpose of one map
`implicit_node.hpp`	`supply_derivatives`, `implicit_value`: a value on the tape carrying rows obtained elsewhere; includes no solver header
`ode_linalg.hpp`, `ode_jacobian.hpp`	dense LU and eigenvalues; the Jacobian from a System’s hook or by forward AD
`ode_substeps.hpp`	`SubsystemSubsteps`: a step that integrates a System’s subsystem (`HasSubsystem`) in substeps of its own, on whichever explicit tableau derives from it
`ode_step_rkck.hpp`, `ode_step_dopri.hpp`, `ode_step_rodas.hpp`	the three steppers; the two explicit ones, `rkck` and `dopri`, record a run, carry `step_adjoint` and substep a subsystem
`ode_callback_system.hpp`	`CallbackSystem`, a System over `std::function`s
`ode_solver_internal.hpp` → `ode_solver.hpp`	`SolverInternal` (stepping, the record), `Solver<System>` (`advance_*`, `recording()`, `sweep_adjoint`, `advance_schedule`)
`sweep.hpp`, `ode_steady_state.hpp`	segment helpers for a consumer driving a sweep itself; Newton to a fixed point with implicit-function-theorem sensitivity
`rcpp_interface_helpers.hpp`, `solver_interface.hpp`	the two headers that need R: `Solver_*_impl` behind the R bindings, including `Solver_loss_and_gradient_impl`

Everything is C++20 (concepts). A consumer compiles with CXX_STD = CXX20 and the two XAD defines from src/Makevars, -DXAD_NO_THREADLOCAL -DXAD_USE_STRONG_INLINE, on every platform: the storage class of the active tape does not change the symbol’s name, so a translation unit built without them reaches the same tape through the other storage class.

Why there is compiled code at all

odelia is almost header-only — the ODE Solver, interpolator, and ODE interface are templates that compile into whichever package includes them. The exception is the XAD automatic-differentiation Tape runtime:

inst/include/XAD/Tape.hpp declares the Tape<T,N> members and, near the bottom, does extern template class Tape<double>; (and friends). extern template tells every translation unit “do not instantiate Tape here — it is compiled elsewhere.”
src/Tape.cpp holds the only copies of the Tape<T,N> member-function definitions (~Tape(), computeAdjoints(), resetTo(), …) and the explicit instantiations (template class Tape<double>;, via the MAKE_TAPE_TLS macros). These compile into odelia.so / odelia.dll and nowhere else.

This split is a deliberate compile-time optimisation (the heavy Tape code is compiled once, in odelia, not in every consumer). Its consequence is that any consumer that instantiates Solver<…> emits undefined references to xad::Tape<double,1>::~Tape() et al. that must be resolved against odelia’s compiled library.

The single-`active_tape_` invariant

Tape.hpp declares static XAD_THREAD_LOCAL Tape* active_tape_;, defined once in Tape.cpp. This is one active-tape pointer per thread, per process — AD recording (getActive(), activate(), deactivate()) all key off it. The whole design assumes a single definition of this symbol in the process.

Keeping Tape compiled in exactly one place (odelia’s library) is what guarantees that single definition. Any change that causes Tape to be instantiated in each consumer (e.g. making it header-only, or shipping the definitions as an includable impl header) gives every consumer DLL its own active_tape_, which is only safe if a given AD computation never crosses the odelia↔︎consumer DLL boundary — a subtle, easy-to-violate condition. Do not do this without understanding the implication.

How downstream packages resolve the symbols, per platform

A package that LinkingTo: odelia and instantiates Solver must get odelia’s compiled Tape symbols resolved. Two things make that work today:

odelia loads its DLL with global symbol visibility in R/zzz.R (.onLoad → odelia_load_dll(local = FALSE); see odelia #26/#28). R runs an imported namespace’s .onLoad before it loads the importing package’s own DLL, so odelia’s symbols are present in the process first.
The consumer must load odelia’s namespace before its own DLL. LinkingTo is C++-only and does not load odelia’s namespace; the consumer therefore needs Imports: odelia and an actual NAMESPACE import so R loads odelia (running the .onLoad above) ahead of the consumer’s useDynLib. plant does this with @importFrom odelia odelia_load_dll (mirroring the @importFrom Rcpp evalCpp idiom). devtools::load_all() happens to load Imports first, but a fresh library() only does so via a real NAMESPACE import.

Platform	Linker behaviour for the undefined `Tape` symbol	Status
macOS	`-undefined dynamic_lookup` (flat namespace): allowed in `.so`, bound lazily at first call; resolved against globally-loaded odelia	✅ works (steps 1–2)
Linux	GNU `ld -shared` allows undefined symbols in `.so`, resolved at `dlopen` against globally-loaded odelia	✅ works (steps 1–2)
Windows	mingw `ld` requires every symbol resolved at link time; a DLL cannot carry undefined imports, so the runtime trick in step 1 cannot help. The consumer must link against odelia’s single compiled DLL at build time	✅ works (consumer `Makevars.win`, below)

On Windows the runtime global-load (step 1) cannot satisfy the linker, so the consumer links plant.dll directly against odelia’s one compiled odelia.dll — which exports the XAD Tape runtime (mingw exports all global symbols by default). This keeps exactly one copy of Tape / active_tape_ in the process (the same single-definition guarantee macOS/Linux get at load time), so it is consistent with the invariant above. The consumer adds a Windows-only src/Makevars.win:

CXX_STD = CXX20
## The two XAD defines must match odelia's src/Makevars on every platform.
PKG_CPPFLAGS = -isystem../inst/include/ -DXAD_NO_THREADLOCAL -DXAD_USE_STRONG_INLINE
## odelia's DLL is at <odelia>/libs/<arch>/odelia.dll. Resolve the base libs dir
## in R, then append R's $(R_ARCH) make-variable (e.g. /x64). Do the arch in make
## -- putting $r_arch inside the Rscript -e quotes lets the shell expand it away.
ODELIA_LIBDIR = $(shell "${R_HOME}/bin/Rscript" -e "cat(system.file('libs', package='odelia'))")
PKG_LIBS = "$(ODELIA_LIBDIR)$(R_ARCH)/odelia.dll"

(plant uses this; verified green on windows-latest R-CMD-check.) odelia must therefore keep exporting these symbols from its DLL — do not add a restrictive .def or -Wl,--exclude-all-symbols to odelia’s build.

A right-hand side handed in at run time

ode_callback_system.hpp holds CallbackSystem, a System whose right-hand side (and optionally Jacobian and validity predicate) are std::functions handed in at construction, so a callable from any language is stepped by the same solver as a compiled system. It is part of the R-free core and is exercised as plain C++ in tests/standalone/; the R adapter in src/r_system.h is a few lines that wrap Rcpp::Functions, and a Python binding would be the same lines over a py::function. Two contracts the adapter relies on and the core keeps: evaluation is lazy (setting the state makes no call; reading the rates does), and the last call of an accepted step is at the accepted state.

The implicit stepper takes its Jacobian from a System’s own ode_jacobian() hook when there is one, else by forward-mode AD on a rebind()-able system (ode_jacobian.hpp); fd_jacobian() there is the one-line finite-difference body for a hook. A system declaring ode_autonomous() is not asked for a df/dt term.

ode_steady_state.hpp reuses that Jacobian and the dense LU away from the stepper: SteadyState<System> solves f(y*) = 0 by damped Newton, reads the fixed point’s stability off the eigenvalues of df/dy (ode_linalg.hpp, Hessenberg reduction then Francis QR), and takes the parameter sensitivity dy*/dθ = −(df/dy)⁻¹ (df/dθ) from the same factorisation. df/dθ is the same forward sweep as df/dy with the direction on a parameter instead of a state component; which parameters is the System’s ad_parameters() hook, a vector of pointers to them in a fixed order, which the rates must read live (a quantity cached from a parameter at construction carries no tangent and yields a zero column; check_parameters() detects that through rebind()). Newton finds any root, attracting or not; solve_with_warmup() integrates the transient when the root it finds is not attracting. Endpoint-only and tape-free by design: the transient is never differentiated, solve() refuses an active scalar type, and a consumer that needs y* on its own adjoint tape attaches the sensitivity rows as a supplied derivative rather than taping the Newton iteration. The callback System reaches the same Newton through its Jacobian hook, which is how ode_steady_state() serves an R right-hand side (src/r_steady_state_interface.cpp); its parameter sensitivity is a finite difference through func done in R, since nothing differentiates R code.

ode_step_dopri.hpp is a third stepper, Dormand-Prince 5(4), whose value over Cash-Karp is a free continuous extension of order 4: SolverInternal::interpolate reads the state anywhere inside the last accepted step from the step’s own stages, and Solver::advance_collect(times, dense = true) uses it to report at requested times without making a step end on each. Under the other steppers the interpolant is cubic Hermite on the step’s endpoints, one order short. Dormand-Prince records, replays and sweeps a run exactly as Cash-Karp does: the same six-entry solved_row, filled with its five later stages and the evaluation at the state the step ends at, and a DopriStep::step_adjoint on its own tableau. The entries mean different stages under the two tableaux, so each step_record names its Method and a replay under the other stepper is refused.

The solver core contains no R

Everything in inst/include/odelia/ except solver_interface.hpp and rcpp_interface_helpers.hpp compiles and runs with no R installation at all. util::stop() throws a std::runtime_error, which Rcpp turns into an ordinary R error at the package boundary; util::warning() writes to std::cerr. This is what lets non-Rcpp consumers use the solver, and what lets consumers that are Rcpp packages — leaf, for one — run their C++ test suites without an R session (traitecoevo/leaf_cpp#11, odelia #43).

tests/standalone/ holds the guard: r_free.cpp, a translation unit that includes the whole core and integrates the Lorenz system (compiled, and through a callback), built by make test-cpp with -I inst/include and nothing else. CI runs it on a runner with no R. It fails loudly if an R header creeps back in, and — the subtler case — if a header comes to rely on a standard-library facility it never includes and was quietly getting from R’s headers. spline.hpp was using assert without <cassert> on exactly that basis.

Beside it, stage_times.cpp is built at -O2 under the compiler’s own contraction default, and fails if a step, a replay and a sweep do not hand a System the same stage times, each rounded as util::time_in_step rounds it.

Note that the standalone build compiles src/Tape.cpp alongside the test. That is the single-Tape-object rule above, not an R dependency: any consumer that instantiates Solver needs the tape runtime from somewhere.

Contract for `LinkingTo: odelia` consumers

Add odelia to LinkingTo: and Imports: in DESCRIPTION.
Compile as C++20 with -DXAD_NO_THREADLOCAL -DXAD_USE_STRONG_INLINE in src/Makevars (and Makevars.win), matching odelia’s own.
Add a real NAMESPACE import from odelia (e.g. @importFrom odelia odelia_load_dll) so odelia’s namespace — and its global-load .onLoad — runs before your package’s useDynLib.
macOS / Linux: nothing else is required (load-time resolution via the global .onLoad above).
Windows: add a src/Makevars.win linking plant.dll against odelia’s DLL (see the snippet in the platform table above).

Do NOT change without coordinating downstream

The extern template / Tape.cpp split, or making Tape header-only — see the single-active_tape_ invariant above.
The set of explicit Tape instantiations in src/Tape.cpp (consumers rely on the exact instantiated types).
XAD_REDUCED_MEMORY in inst/include/XAD/Config.hpp. It sets how the tape stores a recorded operation, so odelia’s compiled Tape and every consumer must be built with the same setting; it lives in the header rather than a -D flag so that a consumer cannot leave it out.
The global-visibility load in R/zzz.R (odelia_load_dll(local = FALSE)) — needed for macOS/Linux.
The set of symbols exported from odelia.dll on Windows — do not add a restrictive .def file or -Wl,--exclude-all-symbols, or Windows consumers can no longer link the Tape runtime.
Reaching for R from anywhere in the core headers (see the section above). Consumers build against them without R.

Any of these will break downstream linking/loading silently (odelia’s own checks will still pass). Coordinate with consumers (plant) and update this document.