Metadata-Version: 2.4
Name: aadc-ode
Version: 0.6.3
Summary: ODE integration with exact AAD gradients on aadc-ng: BDF and RK solvers, kernel-tape replay and off-tape LU solve as compute blocks
License-Expression: LicenseRef-Matlogica-EULA
Classifier: Development Status :: 4 - Beta
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Programming Language :: C++
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS :: MacOS X
Classifier: Operating System :: Microsoft :: Windows
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: Intended Audience :: Science/Research
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE.txt
License-File: THIRD-PARTY-NOTICES.txt
Requires-Dist: aadc<3,>=2.12.0
Requires-Dist: numpy>=1.21.0
Dynamic: license-file

# aadc-ode

ODE integration with **exact** adjoints on aadc-ng.

Two families of solvers:

**C++ kernel replay** (on tape) — the RHS `f(y, p)` is recorded **once**
into its own tape, and each step records a compute block that replays it.

| Solver | Method | For |
|---|---|---|
| `bdf_solve` | Semi-implicit Euler, frozen Jacobian | Stiff systems |
| `rk4_solve` | Classical Runge-Kutta 4th order | Non-stiff, fixed step |

**Python discrete adjoint** — adaptive integration forward, then a backward
sweep computes exact dJ/dp via vector-Jacobian products from the AAD kernel.

| Solver | Method | For |
|---|---|---|
| `adaptive_rk45` | Dormand-Prince 4(5), adaptive step | Non-stiff, high accuracy |
| `fixed_rk4` | Classical RK4, fixed step | Non-stiff |
| `discrete_adjoint` | Backward sweep (one VJP per stage) | Exact dJ/dp for ERK |

**What "exact" means here**, because it is the one place a knowledgeable user
could be misled: exact for the **discrete map with the recorded step
sequence**. The accepted steps `h_list` are frozen and the adaptive controller
is *not* differentiated — that is the semantics of the cited paper, not an
approximation introduced by this port. It is also why the backward sweep needs
the forward solve's time grid and cannot invent one:

```python
t_list, x_list, h_list, k_stages = aadc_ode.adaptive_rk45(f, x0, (t0, tf))
dJdp = aadc_ode.discrete_adjoint(aad_rhs, x_list, h_list, k_stages,
                                 aadc_ode.DORMAND_PRINCE_45, p,
                                 dJdx_T=dJdx_T,
                                 t_list=t_list)     # <- pass it
```

Pass `t_list` (or `t0=` if that is all you have). Omitting both reconstructs
the stage times from **0.0**, which is correct only for an autonomous RHS or a
solve that started at `t=0` — and wrong, silently, otherwise. It now warns;
`t0=0.0` is how an autonomous RHS says so and silences it.

`AadRhs` records the RHS **once** and replays that one recording at every
stage of every step, so a plain `if x[0] > 0:` inside the RHS freezes its
branch at the recording point and uses it at states arbitrarily far away —
wrong values *and* flipped derivatives. Recording therefore **raises** if the
tape reports any active-to-passive conversion; use `iif` / a masked blend, or
`allow_passive=True` if the branch really is constant over the whole
trajectory.

The discrete adjoint is the method from:

> R. Martins, E. Lakshtanov. *A C++ implementation of the discrete adjoint
> sensitivity analysis method for explicit adaptive Runge-Kutta methods
> enabled by automatic adjoint differentiation and SIMD vectorization.*
> Applied Mathematics and Computation, 2025.
> [arXiv:2410.01911](https://arxiv.org/abs/2410.01911)

```
src/bdf.{hpp,cpp}            C++ library: KernelBlock, LUSolveBlock, drivers
src/bdf_selftest.cpp         native checks, incl. the two Python cannot express
bindings/aadcbdf_pybind.cpp  ONE type_caster — where gradients are won or lost
packaging/aadc_ode/          wheel: C++ extension + Python solvers (rk, adjoint)
tests/                       every assertion is on a DERIVATIVE
ci/build.sh · ci/build.ps1   compile + link (no code generation step)
```

## Build and test

```bash
. ci/pins.env && SDK_VER="$SDK_VERSION" SDK_SRC_VERSION="$SDK_SRC_VERSION" \
  CORE_WHEEL_VERSION="$CORE_WHEEL_VERSION" bash ci/build.sh python3.10 python3.11 python3.12 python3.13 python3.14
pip install dist/aadc_ode-*.whl
pytest tests/ -q
```

A requested interpreter that is not installed is a fatal error (matching
`ci/build.ps1`). For a local run that should just build whatever is present,
set `STRICT=0` to skip missing interpreters instead.

Windows: `powershell -File ci\build.ps1 -Pythons 3.10,3.11,3.12,3.13,3.14`, from
inside an MSVC environment. CI covers linux-x86_64, linux-aarch64,
macos-arm64 and windows-x64.

## Use

```python
import aadc, aadc_ode

lam = [1.0, 10.0, 100.0, 1000.0]

# 1. Record the RHS once, in ordinary Python.
k = aadc_ode.Kernel()
states, _ = k.begin(n_states=4)
k.end([-lam[i] * states[i] for i in range(4)])

# 2. Record the loop onto an outer tape. 100 steps -> 100 blocks, ONE kernel.
f = aadc.Functions()
f.start_recording()
y0   = [aadc.idouble(1.0) for _ in range(4)]
args = [v.mark_as_input() for v in y0]
xT   = aadc_ode.bdf_solve(k, y0, dt=1e-3, n_steps=100, jac_lag=10)
cost = sum((x * x for x in xT[1:]), xT[0] * xT[0])
out  = cost.mark_as_output()
f.stop_recording()

# 3. Replay and differentiate, at any point — not just the recording point.
ws = f.create_workspace()
for a in args:
    ws.set_val(a, 1.0)
ws.forward()
ws.set_diff(out, 1.0)
ws.reverse()
print([ws.diff(a) for a in args])
```

`begin` hands back **copies** of the library-owned input variables, and that
direction is the whole design. An aadc-ng tape slot is keyed by the **address**
of its variable, and this binding marshals `Real` by value — so a
`mark_states(y)` taking your list would mark throwaway copies while you kept
using the originals. Zero gradients, nothing raised. With `begin`, your
expressions sit downstream of the real inputs and the adjoint reaches them.

## Design notes

Four properties that affect how you can use this package, rather than how it
came to be:

**Lane-generic blocks.** The compute blocks are `ComputeBlockT`, so a tape
containing them runs on a vector workspace as well as a scalar one. The
`avx2_lanes` self test drives four AVX2 lanes with four different initial
conditions and requires them to reproduce four scalar runs **bitwise**:

```
PASS avx2_lanes    4 lanes vs 4 scalar runs, worst |delta| = 0.000e+00
```

**One reverse sweep per step, not one per state.** A kernel block seeds every
output with the outer adjoint and reverses once, computing `x̄ += Jᵀȳ` directly
instead of assembling a full Jacobian. On a 27-state model that is one sweep
per step rather than 27.

**Thread-safe replay.** Each thread gets its own kernel replay workspace. The
`concurrent_replay` self test runs eight threads against one tape and requires
bitwise agreement with the sequential run.

**The frozen Jacobian is an option, and it is an approximation.** The
semi-implicit step's cached `1 - dt·J_ii` is a **record-time constant**:
replaying with bumped inputs reuses the Jacobian captured while recording.
That is the standard semi-implicit trade — the Jacobian is a preconditioner,
not part of the answer — but it *is* an approximation to the derivative of the
scheme. `freeze_jacobian=False` turns it off; `jac_lag=1` refreshes it every
step.

## What is not differentiated

`luSolve` takes the LU factorisation and pivot vector as **passive** data and
differentiates only the right-hand side. That is the right split for a BDF
step — the iteration matrix is a preconditioner rebuilt from the cached
Jacobian, not a differentiable input — but it is a real limitation for anyone
wanting `dx/dA`. Adding it means differentiating the factorisation itself
(`Ā = -x̄ xᵀ` and friends), which is a separate piece of work.

## Testing philosophy

Every test asserts a **derivative**. A binding that passivates `Real`, a second
libaadc-ng image, or a block whose reverse never reaches the outer tape all
produce values correct to the last digit and gradients of exactly zero — so a
value check proves nothing.

The stiff adjoints are checked against a **closed form**, not finite
differences. By T = 0.1 the λ=1000 component has decayed to ~1e-60, and a
central difference of a quantity that small is pure cancellation noise; FD gets
that adjoint wrong by 0.15%, which would either fail a correct implementation
or force a tolerance loose enough to accept a wrong one.

The `avx2_lanes` and `concurrent_replay` checks are native C++ and bound as
`aadc_ode.selftest()`. They have to be: `aadc_ng`'s own Python
external-function block is scalar-only and its callbacks are
Jacobian-supplying, so a Python-only suite could not exercise either — and
those are precisely the capabilities this package adds.

## Packaging

Per-interpreter (pybind11 is not abi3), on top of the **abi3** `aadc_ng`,
bridged by the `_aadc_ng._C_API` capsule. This wheel ships **no**
libaadc-ng: it links the copy inside the installed `aadc_ng` wheel and finds
it again by RPATH (`$ORIGIN/../aadc_ng`, `@loader_path/../aadc_ng`, and on
Windows by import order — `delvewheel` is deliberately not run).

Recording state is thread-local **inside** libaadc-ng, so a second image gives
the extension and the `aadc` module private recording contexts. For this
package that failure would be total: every entry point either opens a
recording or records a block onto one someone else opened.
`test_exactly_one_library_image` asserts it directly.

Unlike `aadc-example-pybind11` this build does **not** define
`AADCNG_ALLOW_TO_PASSIVE_BOOL`. That macro enables active-to-passive bool
conversions across the whole translation unit, and a package whose entire
purpose is keeping an exact adjoint through a stiff solver should not switch on
implicit passivation to compile itself. `bdf.hpp` includes only
`<aadcNG/idoubleNG.h>` — never the `<aadcNG/aadcNG.h>` umbrella — which keeps
the deleted `iboolNG::operator bool()` doing its job.

## Licence and support

`aadc-ode` is distributed under the MatLogica EULA and depends on the `aadc`
core, which runs in **Community Edition** unless licensed — non-commercial and
academic use only. The full licence and the third-party notices travel inside
the wheel, under `aadc_ode-<version>.dist-info/licenses/`.

For commercial licensing, advanced features or support:
[matlogica.com](https://matlogica.com) · info@matlogica.com
