Metadata-Version: 2.2
Name: yabi-torch
Version: 0.1.0
Summary: A batched, differentiable ODE solver — YABI's BDF + continuous adjoint, wrapped as a torch-native solve_ivp.
Keywords: ode,bdf,torch,differentiable,adjoint,neural-ode
Author: YABI contributors
License: Apache-2.0
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Mathematics
Project-URL: Homepage, https://github.com/dominik-skanda/yabi-bdf
Project-URL: Repository, https://github.com/dominik-skanda/yabi-bdf
Project-URL: Issues, https://github.com/dominik-skanda/yabi-bdf/issues
Requires-Python: >=3.10
Requires-Dist: numpy
Requires-Dist: torch>=2.0
Description-Content-Type: text/markdown

# YABI — Yet Another BDF Integrator

> **Note on origin:** This is a complete rewrite of the BDF integrator developed
> by Dominik Skanda (based on his dissertation work at University of Freiburg).
> The mathematical formulation follows BDF methods with sensitivity generation as presented in the [dissertation by Dominik Skanda](https://freidok.uni-freiburg.de/data/8787).  
> It supports dense (LAPACK), sparse (UMFPACK/KLU via SuiteSparse), and matrix-free Newton–Krylov
> (GMRES/FGMRES/BiCGSTAB via Jacobian–vector products) linear solvers through a pluggable `LinearSolver` interface.
> The rewrite was carried out using **OpenCode** with the **Big Pickle** model and
> **OpenRouter** free models. All algorithmic components (Nordsieck predictor,
> simplified Newton, matrix scaling/conditioning strategy, corrector acceptance,
> and sensitivity replay) were reimplemented to match the mathematical structure
> of the original approach while exposing a clean modern interface.

YABI is a variable-order, variable-step BDF integrator for stiff ordinary
differential equations. It provides forward, adjoint (reverse) and simultaneous
sensitivity analysis through CppAD automatic differentiation.

- **Orders 1–6**, variable order and variable step size.
- **Three right-hand-side forms**: autonomous `f(y)`, time dependent `f(t, y)`,
  and time/parameter dependent `f(t, y, p)`.
- **Sensitivities**: forward mode, adjoint/reverse mode, and simultaneous
  (variational) sensitivities.
- **Stiff problems**: implicit corrector with Newton iteration, Jacobian reuse,
  and LAPACK linear solves.
- **Python bindings**: `import yabi` — zero-copy NumPy I/O with an analytic
  Jacobian (v1, analytic path; CppAD AD and sensitivities in Python are the next
  milestone — see "Python bindings" below).
- **C++17**, BLAS/LAPACK, and an external CppAD installation.

## yabi-torch — `pip install yabi-torch`

A batched, differentiable `solve_ivp` for torch, on top of the BDF + continuous adjoint:

```bash
pip install yabi-torch
```

```python
import torch
from yabi_torch import solve_ivp

def f(y):  # y' = f(y) — any torch.nn.Module or callable
    return torch.stack([-0.04 * y[0] + 1e4 * y[1] * y[2],
                        0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1] ** 2,
                        3e7 * y[1] ** 2])

y0 = torch.tensor([1.0, 0.0, 0.0], dtype=torch.float64, requires_grad=True)
yT = solve_ivp(f, y0, t_span=(0.0, 1e5), rtol=1e-6, atol=1e-10)
yT.sum().backward()   # y0.grad = d(sum y(T)) / d y0 via the continuous adjoint
```

`f` is a `torch.nn.Module` (or any callable `y -> f(y)` on a float64 tensor); the Jacobian is **automatic**
(`torch.autograd`) unless you pass an analytic `jac`. See "Python bindings" below for the full API
(`solver="dense"|"matrix_free"`, `threads=N` for a parallel batch, `device="cuda"`).

## Quick Verification

Prove tolerance adherence, first-order convergence, and higher-order sensitivity
in one command. This regenerates the work–precision diagram
(`work_precision_diagram.png`) and prints the slope summary plus the
higher-order (forward/reverse, orders 1–4) sensitivity output locally.

```sh
cmake -S . -B build && cmake --build build -j
python3 reproduce_benchmarks.py
```

`reproduce_benchmarks.py` runs two self-checking testbench binaries and parses
their output:

A broader, machine-readable performance/accuracy table (runtime, Jacobian
evaluations, factorizations, steps, order, and error vs an independent SciPy
Radau reference) is generated by the telemetry harness:

```sh
cmake --build build --target yabi_telemetry
python3 benchmarks/report.py build/yabi_telemetry
```

This writes `benchmarks/TELEMETRY.md` (the table) and
`benchmarks/work_precision.png` (the error-vs-rtol plot). See `benchmarks/` and
the reference values recorded in `benchmarks/telemetry.cpp`.

- `build/yabi_test_11_tolerance_and_order` — the work–precision proof (tolerance
  adherence + first-order step convergence).
- `build/yabi_test_28_second_order_sensitivity` — exact forward and reverse
  sensitivities of orders 1–4.

Plotting requires `matplotlib` (the build and the two binaries themselves do
not):

```sh
pip install matplotlib
```

## Dependencies

| Dependency | Notes |
|---|---|
| CppAD | Optional. For automatic differentiation / examples; can be disabled or not required in some builds. Must provide `ADFun::capacity_order` (2013 or newer) if used. |
| BLAS | Fortran `dgemm_`, `dgemv_`, `dger_`, `daxpy_`, `dscal_`, `dnrm2_`. |
| LAPACK | Fortran `dgetrf_`, `dgetrs_` (for dense linear solver). |
| SuiteSparse (UMFPACK/KLU) | Optional, for sparse linear solvers via the pluggable `LinearSolver` interface. |
| C++17 compiler | GCC or Clang. |
| CMake ≥ 3.16 | Build only. |
| pybind11 ≥ 3.0 | Optional, for the `yabi` Python module (`-DYABI_BUILD_PYTHON=ON`). `pip install pybind11`. |
| python3-dev | Optional, for the Python module (Python.h headers). |
| matplotlib | Optional. Only for `reproduce_benchmarks.py`'s work–precision plot (and `pyplot/plot_all.py`). |

CppAD is distributed under EPL-2.0 with a GPL-2.0-or-later secondary license.
BLAS and LAPACK are BSD-style. YABI itself is Apache-2.0.

On Debian/Ubuntu:

```sh
sudo apt install build-essential liblapack-dev libblas-dev cmake libcppad-dev
# optional, for the Python bindings (T37):
sudo apt install python3-dev && pip install pybind11
```

## Layout

```
yabi-bdf/
  include/yabi/
    bdf_integrator.hpp        public class declaration
    linear_algebra.hpp        BLAS/LAPACK bindings
    detail/
      core.ipp                construction, reset, configuration, results
      step_control.ipp        step size, predictor, order control, output
      newton.ipp              Newton iteration and corrector acceptance
      sensitivity_setup.ipp   sensitivity storage and tape lifetime
      sensitivity_forward.ipp forward-mode replay
      sensitivity_reverse.ipp reverse/adjoint replay
  src/
    bdf_integrator.cpp        explicit instantiation for double
  examples/
  tests/
    unit/                     isolated component unit tests
    component/                cross-component integration tests
    system/                   mathematical benchmarks and regressions
```

The `.ipp` files contain the member definitions of the class template. They
are implementation details: include `<yabi/bdf_integrator.hpp>`, not the
`.ipp` files.

## Build

```sh
cmake -S . -B build [-DCPPAD_ROOT=/path/to/cppad]
cmake --build build -j
ctest --test-dir build --output-on-failure
```

If CppAD is not installed in standard paths:

```sh
cmake -S . -B build -DCPPAD_ROOT=/path/to/cppad
cmake --build build -j$(nproc)
ctest --test-dir build
```

When CppAD is installed system-wide, `-DCPPAD_ROOT` / `CPPAD_ROOT` can be
omitted. CMake first tries `pkg-config`, then looks for `cppad/cppad.hpp` and
`libcppad_lib`.

## Python bindings (T37)

An optional `yabi` Python module exposes the integrator on the analytic path:
the RHS and analytic Jacobian are Python callables, and the state buffers are
zero-copy NumPy arrays.

```sh
pip install pybind11         # >= 3.0 (Python 3.14)
sudo apt install python3-dev # Python.h headers
cmake -S . -B build-py -DYABI_BUILD_PYTHON=ON
cmake --build build-py --target yabi_python
ctest --test-dir build-py -R python_bindings
```

```python
import numpy as np, yabi
def f(y):   return np.array([...])   # RHS, n doubles
def jac(y): return np.array([...])   # column-major J, n*n doubles
it = yabi.BDFIntegrator(f, jac, n)
y0 = np.array([...]); atol = np.array([...])
it.reset_integration(y0, max_step, rtol, atol)   # zero-copy y0
it.set_algorithm_settings(5, np.finfo(float).eps, True, 1)
it.integrate(t_end)
it.get_solution(y_out)                            # zero-copy write into y_out
```

`bindings/python/test_yabi.py` is the Brusselator smoke test;
`bindings/python/test_autograd.py` checks the continuous-adjoint binding against
central differences.

The first-order **continuous adjoint** (stable reverse sensitivity) is exposed as
`integrator.continuous_adjoint(seed_weights, adj_y0)` — the backward hook for
autograd wrappers.  `bindings/python/torch_wrapper.py` wraps it in a
`torch.autograd.Function` (forward solve + adjoint backward), so gradients flow
through the ODE solve w.r.t. the initial condition.  `bindings/python/yabi_torch.py`
provides a batched, differentiable `solve_ivp(f, y0, t_span, rtol, atol, jac=None,
solver="dense", threads=1, device="cpu")` on top of it: `f` is a `torch.nn.Module`
(or any torch callable `y -> f(y)`), the Jacobian is **automatic** from
`torch.autograd` unless an analytic `jac` is given, and `solver="matrix_free"`
drives the Newton solve with a `torch.autograd` JVP (forward-mode, one `J·v` per
Krylov iteration) instead of the full n×n Jacobian.  `y0` may be a batch; the
gradient is the continuous adjoint.  `threads=N` solves the batch's independent
trajectories concurrently (the binding releases the GIL during the C++ solve and
re-acquires it only for the Python callbacks), bit-identical to serial.
`device="cuda"` offloads the dense Newton LU to cuSOLVER (`CudaLinearSolver`) and,
when the GPU has compute capability ≥ 7.5, evaluates `f`/the Jacobian on-device
(torch CUDA); it requires a CUDA build of the module (`YABI_USE_CUDA=ON`).  Honest
scope: the BDF state stays host-side (host↔device copy per evaluation); a full
device-resident-state BDF is a larger follow-up.  The cuSOLVER LU offload is
verified bit-identical to CPU; the f-on-GPU path is implemented but unverified on
the development GPU (a Pascal GTX 1050, CC 6.1, which torch's cu12x wheels no
longer support).
`test_solve_ivp.py` runs the `torch.autograd.gradcheck`, batch (forward + adjoint),
and stiff Robertson checks (dense + matrix-free); `test_nn_module.py` checks an
`nn.Module` `f` against the matrix exponential and cross-mode agreement;
`test_batch_parallel.py` checks `threads=N` is bit-identical to serial;
`test_cuda.py` checks the cuSOLVER LU offload and (gated on CC ≥ 7.5) the f-on-GPU
path.  Not yet exposed in Python: forward-mode sensitivities, time-dependent RHS,
and the CppAD auto-diff path.

## Usage

### User ODE formulation
- **CppAD / automatic differentiation:** Pass an `ADFun<double>` (tape) describing the RHS — sensitivities/Jacobians can be handled via AD/Internal Numerical Differentiation.
- **Analytical derivatives:** Use the user-ODE interfaces in `include/yabi/user_ode.hpp` and the manual sensitivity examples (`examples/manual_sensitivity_demo.cpp`, `examples/manual_sensitivity_pure.cpp`) to provide analytic Jacobians/sensitivity RHS without requiring CppAD.
- **Custom sparse solvers:** Implement `yabi::LinearSolver` (`include/yabi/linear_solver.hpp`) to plug your own sparse direct/iterative solver; set it via `integrator.setLinearSolver(...)`. UMFPACK/KLU are provided via SuiteSparse; dense LAPACK is the default.

```cpp
#include <yabi/bdf_integrator.hpp>

CppAD::vector<CppAD::AD<double>> x(1);
x[0] = 1.0;                       // y(0)
CppAD::Independent(x);

CppAD::vector<CppAD::AD<double>> f(1);
f[0] = -0.5 * x[0];

CppAD::ADFun<double> rhs;
rhs.Dependent(x, f);
rhs.optimize();

yabi::BDFIntegrator<double> integrator(rhs);

double y0[1]       = {1.0};
double abs_error[1] = {1e-12};
double y[1];

integrator.resetIntegration(y0, 0.1, 1e-10, abs_error,
                            false, false, nullptr, nullptr, nullptr, nullptr, 0);
integrator.setAlgorithmSettings(5, std::numeric_limits<double>::epsilon(), true, 1);
integrator.performIntegration(2.0, 100000);
integrator.getSolutionVector(y);
```

Dense output: `getSolutionAt(t, y)` evaluates the continuous extension at any `t`
inside the last accepted step `[t_end - h, t_end]` (it returns `false` outside):

```cpp
double h = integrator.getCurrentStepSize();   // last step size
integrator.getSolutionAt(2.0 - 0.5 * h, y);   // interpolate inside the last step
```

The seven trailing sensitivity-seed arguments can be replaced by a
`yabi::IntegrationOptions` struct (defaults reproduce the `false, false,
nullptr, …, 0` tail):

```cpp
yabi::IntegrationOptions<double> opts;             // taping=false, direct=false,
integrator.resetIntegration(y0, 0.1, 1e-10, abs_error, opts);  // seeds=nullptr, numOfderivatives=0
```

For sensitivity, fill the fields and pass the same struct to any of the three
overloads (autonomous / time / time+parameter):

```cpp
opts.taping = false;
opts.direct = false;
opts.y_seed_list              = yseed_ptr.data();   // double** (per-direction state seeds)
opts.p_seed_list              = pseed_ptr.data();   // double** (per-direction parameter seeds)
opts.relative_error_derivatives = rel_err.data();   // double*  (per-direction rel. tol)
opts.absolut_error_derivatives   = abs_ptr.data();  // double** (per-direction abs. tol)
opts.numOfderivatives         = ndir;
integrator.resetIntegration(y0, t0, p, max_step, 1e-10, abs_error, opts);
```

The RHS type is inferred from the dimensions of the CppAD tape:

| Tape domain | Tape range | Meaning |
|---|---|---|
| `n` | `n` | autonomous: `f(y)` |
| `n + 1` | `n` | time dependent: `f(t, y)` |
| `n + 1 + p` | `n` | time and parameter dependent: `f(t, y, p)` |

`setAlgorithmSettings(max_order, machine_precision, hit_end_time,
error_weight_index)` controls the maximum BDF order (1–6), the precision used
for the final-step test, whether the integrator must land exactly on the end
time, and which error-weight norm to use. All four arguments are validated
before any state is changed.

`max_order` is a **cap**, not a forced order: the controller chooses the order
actually used. Use `getCurrentOrder()` to read it. Likewise `max_StepLength` caps
the step but is not an accuracy knob. Control accuracy with `rtol`/`atol`; because the
controller changes order as the step cap changes, the error is not monotone in
`max_StepLength` (measured examples in `tests/system/README.md`).

To **force** a fixed order (e.g. for per-order convergence measurement), use
`setFixedOrder(q)`:

```cpp
integrator.setFixedOrder(4);        // pin the order to 4 (ramps up 1..4, then stays)
integrator.getCurrentOrder();       // read the order actually used
```

### Step-controller adaptation (stability)

The step controller is the original 2011 BDF algorithm with one small, deliberate
addition for stability: on a rejected step it estimates the local truncation error at the
next-lower order (from the highest Nordsieck term) and reduces the order when that gives a
smaller or already-acceptable error. This order reduction lets the
integrator recover from stiff-oscillatory problems — e.g. the Oregonator's relaxation
oscillations, where high-order BDF is unstable at *any* step size. It is a **no-op for
well-behaved problems** (they accept on the first attempt, so no order reduction ever fires),
and each reduction is recorded in the sensitivity tape so forward/reverse replays reconstruct
the correct order sequence.

### Extreme-stiffness behaviour (EMEP)

On the hardest stiff benchmarks (the CWI "Test Set for IVP Solvers" **EMEP** problem — 66
species, stiffness ratio ~1e14) the plain defaults used to stall with `too much steps rejected`.
The cause was a one-character boundary bug in the step controller.

The minimal-step floor sets `h = 4·eps·|t|` *exactly* (`eps` = machine precision), yet the
"accept the step at minimal size" escape tested `h < 4·eps·|t|`.  The strict inequality never
matched the floored value, so once a stiff transient drove the step down to the floor, the
local-error test could no longer be satisfied and the step was rejected forever.  The fix tests
`h <= 4·eps·|t|` (and `h > 4·eps·|t|` in the still-reducible branch), so at the floor the step
is accepted with a warning rather than rejected indefinitely — the standard `hmin` handling
(Hairer & Wanner §III.7).

With that, EMEP matches its CWI reference to ~8e-6 with **default Newton settings and no fixed
order** — only the standard stiff-BDF max order of 5 (BDF6 is not A(α)-stable; see
`tests/system/45_emep.cpp`).

`setNewtonLineSearch(true)` remains an independent, opt-in feature: a **damped Newton**
(backtracking line search) that halves the Newton correction while the full step would increase
the BDF residual — the standard fix for Newton overshoot (Nocedal & Wright §11; Hairer & Wanner
§IV.8).  It is a **mathematical no-op where Newton does not overshoot** (the full step is
accepted), so it is **off by default** — the one extra residual evaluation reorders
floating-point, so the BASELINE anchors stay byte-identical only with it off.  EMEP does not
need it; it is kept for general Newton-overshoot robustness.

References: Hairer, Nørsett & Wanner, *Solving ODEs I/II*; Gear (1971), *Numerical Initial Value
Problems in ODEs*; Byrne & Hindmarsh (1975); Nocedal & Wright, *Numerical Optimization* §11;
Mazzia et al., "Test Set for IVP Solvers".

### Relationship to the dissertation (method provenance)

The mathematical formulation is the [dissertation by Dominik Skanda](https://freidok.uni-freiburg.de/data/8787)
(University of Freiburg).  **That method is correct and remains the baseline**: a variable-order,
variable-step BDF corrector with a Milne (corrector–predictor) local-error estimator and an
**error-per-step (EPS)** acceptance test — the dissertation's original estimator, which is *kept*.

Two **deliberate, literature-referenced robustness additions** sit on top of that baseline.  Both are
**no-ops on well-behaved problems** (they fire only on the pathological cases they exist for), so the
BASELINE anchors are byte-identical to the dissertation method for the standard benchmarks:

1. **Order reduction on a rejected step** — *[§ "Step-controller adaptation" above]*.  Source:
   Gear (1971); Byrne & Hindmarsh (1975).
2. **Damped Newton** (backtracking line search, opt-in) — *[§ "Extreme-stiffness behaviour" above]*.
   Source: Nocedal & Wright, *Numerical Optimization* §11; Hairer & Wanner §IV.8.

A third change is a **bug fix, not an addition**: the "accept the step at minimal size" escape
that was present in the original code (since the initial commit) tested `h < 4·eps·|t|` while the
minimal-step floor sets `h = 4·eps·|t|` *exactly*.  Correcting the test to `h <= 4·eps·|t|` makes
the escape fire as it was always intended, which is what lets the extreme-stiffness EMEP problem
complete (see above).  It is likewise a no-op on well-behaved problems: a healthy step is far above
the floor, so the comparison outcome is unchanged.

## Sensitivity API

### Forward mode

Enable taping in `resetIntegration`, integrate to the output time, then call
`forward`:

```cpp
double y_seed[2] = {1.0, 0.0};    // d/dy0
double p_seed[2] = {0.0, 0.0};
double direction[2];

integrator.forward(1, y_seed, p_seed, false);
integrator.getForwardSolutionVector(direction);
```

For a parameter derivative, set `y_seed` to zero and the corresponding entry
of `p_seed` to one. Calling `forward(..., true)` at the same output time resets
the forward replay so that a new direction can be requested.

### Reverse / adjoint mode

Integrate to the output time with taping enabled, then seed the adjoint with
the weights of the scalar output of interest:

```cpp
double weight[2] = {1.0, 0.0};
double* adj_species[2];
double* adj_parameter[2];

integrator.reverse(1, weight);
integrator.getReverseSolutionVector(adj_species, adj_parameter);
```

`adj_species[i][k]` is the k-th order adjoint of state `i`;
`adj_parameter[i][k]` is the corresponding parameter adjoint.

### Continuous adjoint (first order, stable reverse)

The discrete adjoint `reverse()` above is exact for the discretization but
A-stability-bound over long stiff intervals (see T24). `continuousAdjoint()`
integrates the adjoint ODE `dλ/dt = -J(y)ᵀλ` backward in reversed time with the
same BDF integrator, so it stays stable where the discrete replay diverges
(e.g. HIRES to steady state):

```cpp
double weight[2] = {1.0, 0.0};   // d g / d y(T)  for  g = w^T y(T)
double adj_y0[2];
double adj_p[1];                 // one entry per parameter (RHS_type 2)
integrator.continuousAdjoint(weight, adj_y0, adj_p);
// adj_y0[i] = d g / d y0_i,   adj_p[j] = d g / d p_j
```

It requires only that `performIntegration` has run (the forward trajectory is
recorded automatically) and works with both CppAD (AD `Jᵀλ`/`(∂f/∂p)ᵀλ` via
`Reverse`) and the pure `UserODERHS` backend (analytic Jacobian). The parameter
sensitivity `∂g/∂p` is accumulated with an augmented backward state (the `μ`
accumulator), and the pure backend supports time- and parameter-dependent RHS
via `setTimeDerivative`/`setParameterDerivative`.

The **second-order** continuous adjoint is `secondOrderAdjoint` — it returns the
Hessian-vector product `H g · s` of `g = wᵀ y(t_end)` in a seed direction `s`,
stably for long/stiff problems (the stable complement to `reverse(2)`):

```cpp
double weights[2] = {1.0, 0.0}, direction[2] = {1.0, 0.0}, hessian[2];
integrator.secondOrderAdjoint(weights, direction, hessian);
// hessian[i] = sum_j d^2 g / d y0_i d y0_j * direction[j]
```

It needs a Hessian-in-direction hook: CppAD provides it via a `Reverse(2)`
sweep automatically; the pure backend adds `setHessianVectorProduct(hv)`.

The **second-order forward** is `secondOrderForward` — it returns the
second-order directional derivative `d²y/dτ²` in a seed direction, stably for
long/stiff problems (the stable complement to `forward(2)`):

```cpp
double direction[2] = {1.0, 0.0}, second[2];
integrator.secondOrderForward(direction, second);
// second[i] = d^2 y_i / d tau^2  along `direction`
```

It needs a second-order directional-derivative hook: CppAD provides it via a
`Forward(2)` sweep automatically; the pure backend adds
`setSecondOrderDirectional(sod)`.

Higher-order forward sensitivity is `forwardOrder(k, seed, result)` — it returns
`d^k y/dτ^k` for arbitrary `k >= 1`, stably for long/stiff problems (the stable
complement to `forward(k)`, which is A-stability-bound at order >= 3):

```cpp
double direction[2] = {1.0, 0.0}, kth[2];
integrator.forwardOrder(3, direction, kth);
// kth[i] = d^3 y_i / d tau^3  along `direction`  (k = 2 coincides with secondOrderForward)
```

For `k >= 2` it needs the k-th-order directional *source* of the RHS via
`AbstractRHS::orderKDirectional`: CppAD provides it automatically (full-Taylor
`Forward`); the pure backend adds `setOrderKDirectional(okd)` and `forwardOrder`
fails loudly without it.  No finite differences anywhere — the pure hook is
exact (analytic), and `tools/generate_derivatives.py` can emit it for you
(see below).

### Generating higher-order derivatives (pure mode, optional tool)

`tools/generate_derivatives.py` is a **decoupled, optional** SymPy code generator
for the exact higher-order sources.  It is not a build dependency and the solver
never calls it — it emits plain C++ in the `UserODERHS` hook signatures, which you
wire in by hand:

```sh
python3 tools/generate_derivatives.py            # print the generated header
python3 tools/generate_derivatives.py --check    # verify the sources numerically
```

```cpp
#include "generated_myode.hpp"   // the emitted header
yabi::UserODERHS<double> rhs(n, generated::f, generated::jacobian);
rhs.setOrderKDirectional(generated::orderKDirectional);
```

Edit the `N` / `Y` / `F` / `MAX_ORDER` block at the top of the script to describe
your ODE, then regenerate.  Requires SymPy (`pip install --user
--break-system-packages sympy`; see PEP 668).  `tests/system/generated_cubic.hpp`
is a committed example of its output, exercised by
`tests/system/39_generated_higher_order.cpp`.

### Simultaneous (variational) sensitivities

Pass seed directions and error controls directly to `resetIntegration`:

```cpp
double* y_seeds[ndir];
double* p_seeds[ndir];
double* abs_seeds[ndir];
double rel_seed = 1e-9;
// fill y_seeds[d], p_seeds[d] and abs_seeds[d] ...

integrator.resetIntegration(y0, t0, p, max_step, rel_error, abs_error,
                            false, false, y_seeds, p_seeds, &rel_seed,
                            abs_seeds, ndir);
integrator.performIntegration(t_end, max_steps);
integrator.getSimultaneousDerivativesVector(solution_derivatives);
```

`examples/sensitivity_demo.cpp` demonstrates all three modes and checks them
against central finite differences.

The forward/simultaneous sensitivity directions are independent, so their Newton
solve + update parallelize across threads (T35):

```cpp
integrator.setSensitivityThreads(4); // 1 = serial (default); >1 = parallel
```

The direction *residual* (CppAD `Forward(1)` tape replay) stays serial — the tape
is stateful — and the parallel solve requires a thread-safe linear solver (the
default dense LAPACK solve is). `tests/component/parallel_sensitivity.cpp` and
`tests/system/48_parallel_sensitivity.cpp` verify parallel == serial bit-for-bit.

### Higher-order (forward / reverse)

`forward` and `reverse` take an `order` argument and support arbitrary order
(1, 2, 3, …). Order 1 uses the legacy single-coefficient replay. Orders 2 and
up use a full Taylor-polynomial replay of the frozen BDF algorithm: every
intermediate quantity is a truncated Taylor polynomial in the sensitivity
parameter τ, the nonlinear RHS is evaluated with CppAD's full-Taylor
`Forward(q, xq)` (`xq.size() == Domain()*(q+1)`), and the returned value is the
order-th directional derivative `d^order y / d seed^order` (i.e. `order!` times
the order-th Taylor coefficient — the *pure* derivative, not the raw Taylor
coefficient).

```cpp
// A single call computes order q directly; no forward(1), forward(2), ...
// warm-up sequence is needed for q >= 2.
integrator.forward(q, y_seed, p_seed, true);
integrator.getForwardSolutionVector(dir);    // d^q y / d seed^q

integrator.reverse(q, weights);              // q-th-order adjoint
integrator.getReverseSolutionVector(adj_sp, adj_p);   // adj_sp[i][k], k = 0..q-1
```

`reverse(q)` returns the adjoint coefficients of orders 1..q: `adj_sp[i][k]` is
the k-th-order directional sensitivity of `sum_j weights[j] * y_j` with respect
to initial state `i` (and likewise for `adj_p`). It is aligned to the same
pure-derivative convention as `forward` (each order `k` is rescaled by `k!` to
undo CppAD's `1/k!` Taylor weighting). The reverse is the *gradient-level* dual
of the forward solution-level derivative: for a scalar ODE with unit seed and
unit weight, `forward(q)` and `reverse(q)[q-1]` coincide.

Limitations:

- `forward(order >= 2)` requires an RHS backend with full-Taylor support
  (the CppAD adapter). A `UserODERHS`-style backend exposes only
  single-coefficient `Forward`, so order ≥ 2 fails loudly with
  `YABI_USE_CPPAD=OFF`.
- `direct` sensitivity mode has no full-Taylor state of the implicit solve and
  stays limited to order 1; higher orders fail loudly.
- `reverse(q)` must follow the `forward(q)` that produced the output state
  (it replays the forward tangent direction), as in the examples above.

The `tests/system/28_second_order_sensitivity.cpp` driver verifies orders 1–4
against the exact solution of `y' = y²`, and
`tests/system/33_higher_order_sensitivity.cpp` covers a coupled nonlinear
system, a non-autonomous system, and a parameter-dependent system.

### Choosing a sensitivity method

YABI exposes **four pure sensitivity methods**; no automatic switching or
method mixing happens — the user selects exactly one. They split along two
axes: *forward vs reverse* (what the seed/weight means) and *discrete vs
continuous* (how the sensitivity is propagated).

| Method | Direction | Propagation | Stable on long/stiff? |
|---|---|---|---|
| `forward(q)` | forward (state) | discrete IND (replay of frozen step Jacobians) | order 1 yes; order ≥ 3 diverges |
| `reverse(q)` | reverse (adjoint) | discrete IND | order 1 yes; order ≥ 3 diverges |
| `continuousAdjoint` / `secondOrderAdjoint` | reverse | continuous ODE `λ' = -Jᵀλ` (re-integrated) | **yes, orders 1–2** |
| simultaneous seeds / `secondOrderForward` / `forwardOrder(k)` | forward | continuous ODE `S' = J·S` (+ higher-order source) | **yes, arbitrary order** |

The discrete (IND) methods are *exact for the discretization* and are the
fastest and most accurate choice on short intervals — but they replay frozen
BDF linearizations, so at order 3–6 over long, stiff intervals the replay is
A-stability-bound and can diverge to huge finite values (T24). The continuous
methods reformulate the sensitivity as another ODE and hand it to the same
stable BDF, inheriting its stability for free; they cost a second (or third)
full integration.

**Decision table**

| Your situation | Use |
|---|---|
| General, short interval, need accuracy + speed | `forward(1)` / `reverse(1)` (IND) |
| Long/stiff interval, forward sensitivity | simultaneous seeds (1st) or `secondOrderForward`/`forwardOrder(k)` (2nd+) |
| Long/stiff interval, reverse sensitivity | `continuousAdjoint` (1st) or `secondOrderAdjoint` (2nd) |
| Higher-order (≥ 3) derivative on a short interval | `forward(q)` / `reverse(q)` (full-Taylor replay) |
| Higher-order (≥ 3) derivative, long/stiff, forward | `forwardOrder(k)` (continuous, stable) |

**Backend-agnostic contract** — every method works identically under CppAD
(`YABI_USE_CPPAD=ON`) and pure mode (`YABI_USE_CPPAD=OFF`):

- 1st-order continuous/simultaneous: CppAD provides `J` (and `∂f/∂p`) via AD;
  pure `UserODERHS` supplies an analytic Jacobian (plus
  `setTimeDerivative`/`setParameterDerivative` for non-autonomous RHS).
- 2nd-order continuous/simultaneous: CppAD provides `H·v` via a `Reverse(2)`
  or `Forward(2)` sweep; pure mode supplies the analytic
  `setHessianVectorProduct(hv)` (adjoint) or `setSecondOrderDirectional(sod)`
  (forward) hooks.
- Higher-order (≥ 3) forward: CppAD provides the k-th-order source via the
  full-Taylor `Forward` automatically; pure mode supplies the analytic
  `setOrderKDirectional(okd)` hook (or generates it — see above).
- No finite differences are used anywhere. The deliberate gaps are: higher-order
  *continuous adjoint* (≥ 3) is not provided (nothing consumes a 3rd-order
  adjoint), and the *discrete* higher-order replay (`forward(q)`/`reverse(q)`,
  `q ≥ 2`) requires the CppAD full-Taylor backend and fails loudly in pure mode.

**Rejected alternatives** (recorded so they are not re-introduced): no
auto-patching of IND (a diverged replay is reported, not silently corrected),
no on-the-fly conservative-mode switching, and no finite-difference fallback.

## Tests

The tests are split into three ctest tiers below `tests/`:

- `tests/unit/` — one focused executable per component (`IntegratorConfig`,
  `ReplayTape`, `NordsieckState`, `StepController`, `NewtonSolver`,
  `SensitivityEngine`, `AdjointEngine`, `DenseLUSolver`/CSC reference path,
  `SparsityPattern`, `CscMatrix`, CppAD RHS sparsity/`J·v`/batch `J·V`,
  `UserODERHS` analytic sparsity/`J·v`/batch `J·V`, the Krylov solver, the
  pluggable-preconditioner/forcing contracts, and the sensitivity-linearization
  seams). These are fast, dependency-light contract tests.
- `tests/component/` — cross-component integration: native CSC assembly and
  the CppAD sparse Jacobian path (`native_sparse`), matrix-free Newton–Krylov
  against the direct HIRES trajectory (`matrix_free_newton`), the same
  matrix-free path through the pure `UserODERHS` backend
  (`matrix_free_pure`, both configurations), a fully J-free solve with no
  sparsity pattern and no CSC assembly (`matrix_free_jfree`, both
  configurations), plus symbolic reuse with UMFPACK/KLU
  (`suitesparse_sparse`, registered only when SuiteSparse is available).
- `tests/system/` — the mathematical benchmark suite: linear problems,
  Robertson, Van der Pol, Brusselator, HIRES, Oregonator, Prothero–Robinson,
  damped pendulum, E5, large sparse/network systems, forward/reverse/
  simultaneous sensitivities, higher-order (forward/reverse, orders 1–4),
  dense output, fixed order, and failure/tolerance regressions.

Every test checks its result against an independent reference solution.
Run a single tier with `ctest --test-dir build -L unit`, `-L component`, or
`-L system`; `ctest` also registers the three self-checking examples under the
`example` label (`example_basic`, `example_nonautonomous`,
`example_sensitivity_demo`).

See `tests/system/README.md` for the measured convergence slopes; the per-order
convergence test for orders 2–6 is now expressible via `setFixedOrder` +
`getCurrentOrder` (`tests/system/32_fixed_order.cpp`).

## License

Apache License 2.0. See [LICENSE](LICENSE) and [NOTICE](NOTICE).


## API Quick Start

### Basic integration (autonomous)
```cpp
#include <yabi/bdf_integrator.hpp>
#include <cppad/cppad.hpp>

int main() {
    CppAD::ADFun<double> rhs;
    CppAD::vector<CppAD::AD<double>> x(1);
    x[0] = 1.0;
    CppAD::Independent(x);
    CppAD::vector<CppAD::AD<double>> f(1);
    f[0] = -x[0];
    rhs.Dependent(x, f);
    rhs.optimize();

    yabi::BDFIntegrator<double> integrator(rhs);
    double y0[1] = {1.0}, abs_err[1] = {1e-10};
    integrator.resetIntegration(y0, 2.0, 1e-10, abs_err, false, false, NULL,NULL,NULL,NULL,0);
    integrator.setAlgorithmSettings(5, 1e-16, true, 1);
    integrator.performIntegration(2.0, 100000);
    double y[1]; integrator.getSolutionVector(y);
    return 0;
}
```

### Forward/Reverse/Simultaneous sensitivities
See `examples/sensitivity_demo.cpp` for full usage of `forward()`, `reverse()`, and simultaneous seeds via `resetIntegration` with seed arrays.

### Pluggable linear solvers (sparse)
```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/detail/umfpack_linear_solver.ipp>

yabi::BDFIntegrator<double> integrator(rhs);
integrator.setLinearSolver(new yabi::UmfpackLinearSolver<double>()); // takes ownership
```
Defaults to dense LAPACK. Supports UMFPACK/KLU via SuiteSparse. See `examples/heat1d_umfpack_literature.cpp`.

### GPU linear solver (experimental, T36/T36b/T36c)
An opt-in cuSOLVER dense-LU solver offloads the Newton `factorize`/`solve` to a
GPU through the same `setLinearSolver` seam:

```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/cuda_linear_solver.hpp>

integrator.setLinearSolver(new yabi::CudaLinearSolver<double>());
```

Build it with `-DYABI_USE_CUDA=ON` (default OFF — the regular build stays
CUDA-free, with zero CUDA headers/libraries). It uses only the host-side
cuBLAS/cuSOLVER C API (no `nvcc`), and links `-lcusolver -lcudart`.

Each instance owns its own CUDA stream + cuSOLVER handle (`cusolverDnSetStream`)
plus pinned host staging buffers (`cudaMallocHost`) with stream-ordered
`cudaMemcpyAsync` copies, so transfer can overlap compute and **multiple
instances run in parallel without serializing** on the default stream
(`tests/component/cuda_parallel.cpp` runs 4 concurrent integrators, each matching
the CPU reference). `tests/component/cuda_solver.cpp` cross-checks the GPU state
solve against dense LAPACK on HIRES to ~3e-16.

**Scope (honest):** this offloads the *linear solve* only — the Jacobian `J` is
still assembled on the host. The sparse paths are the cuDSS adapter (T36d) and the
GPU Krylov solver (T36e, below). Throughput on a small card (GTX 1050) is modest:
~1.2–1.3× for 4 concurrent instances at n ≈ 1000–2000, slower at n ≈ 400
(contention). The cuBLAS/cuSOLVER API is stable across CUDA 12.6 → 13, so the
solver is source-identical on both (only `-arch`/driver change). Requires CUDA
(built against 12.6) + an NVIDIA driver.

### GPU sparse iterative Krylov solver (own implementation, T36e)

For large sparse Newton systems (n ≫ 100) where the dense cuSOLVER factorization
is too costly, `CudaKrylovSolver` runs a BiCGStab Krylov solve on the GPU — cuBLAS
for the vector ops and cuSPARSE SpMV for `P·v` — behind the same `setLinearSolver`
seam. It is **our own implementation** over CUDA primitives (no external library):

```cpp
#include <yabi/cuda_krylov_solver.hpp>
auto *solver = new yabi::CudaKrylovSolver<double>();
solver->setTolerance(1e-10, 0.0, 200, 30);
solver->setForcingTerm(yabi::eisenstatWalkerForcing<double>()); // inexact Newton
integrator.setLinearSolver(solver);
```

Built with `-DYABI_USE_CUDA=ON` (links `-lcusparse -lcublas -lcudart`, no nvcc;
the default build stays CUDA-free). `tests/component/cuda_krylov_solver.cpp`
cross-checks it against dense LAPACK on HIRES to ~3e-16;
`tests/system/cuda_krylov_heat.cpp` solves the method-of-lines heat equation
(n=500) to the analytic solution.

### External sparse solvers (open architecture, T36d)
The solver architecture is **open**: `setLinearSolver` takes any `LinearSolver`,
and `include/yabi/external_sparse_solver.hpp` provides `ExternalSparseSolver` —
the reusable, solver-agnostic half (CSC→CSR handoff + factorize/solve plumbing).
A concrete adapter implements only `buildCsc` (build the solver from CSR) and
`solveOne` (solve one RHS), so adding a solver is a ~30-line adapter that never
touches the integrator core. We **integrate sophisticated open-source solvers,
we do not reimplement them**.

The first adapter is **AMGCL** (MIT, header-only AMG + Krylov):

```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/amgcl_linear_solver.hpp>

integrator.setLinearSolver(new yabi::AmgclLinearSolver<double>()); // BiCGStab + ILU(k)
```

Build with `-DYABI_USE_AMGCL=ON` (default OFF — dependency-free). Provision the
headers once (no sudo): AMGCL → `~/.local/opt/amgcl`, Boost → `~/.local/opt/boost`
(download the AMGCL tarball and `libboost1.83-dev` via `apt download`, extract).
`tests/component/amgcl_solver.cpp` cross-checks it vs dense LAPACK on HIRES to
`4.4e-16`. The BDF Newton matrix is non-symmetric, so the default is
BiCGStab+ILU; AMG+CG is the SPD-only alternative.

The second adapter is **Ginkgo** (BSD, many-core OpenMP/CUDA/HIP):

```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/ginkgo_linear_solver.hpp>

integrator.setLinearSolver(new yabi::GinkgoLinearSolver<double>()); // BiCGStab (OpenMP)
```

Build with `-DYABI_USE_GINKGO=ON`. Provision Ginkgo from its official source
(cmake + install to `~/.local/opt/ginkgo`, OpenMP backend). It uses the OpenMP
(CPU) backend, so it works without a GPU; the CUDA backend is the same adapter
with a `CudaExecutor` (and a supported GPU arch). `tests/component/ginkgo_solver.cpp`
cross-checks it to `1.1e-16`.

The third adapter is **cuDSS** (NVIDIA's sparse *direct* solver, via the
`nvidia-cudss-cu12` wheel):

```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/cudss_linear_solver.hpp>

integrator.setLinearSolver(new yabi::CudssLinearSolver<double>()); // general sparse LU
```

Build with `-DYABI_USE_CUDSS=ON`. Provision it with `pip download
nvidia-cudss-cu12` + extract to `~/.local/opt/cudss` (never vendored). cuDSS has
real `ANALYSIS/FACTORIZATION/SOLVE` phases, so it matches the BDF
factorize-once/solve-many pattern; `tests/component/cudss_solver.cpp`
cross-checks it to `1.7e-18`. All three adapters reuse the same
`ExternalSparseSolver` base — the libraries are **never vendored**, only the
thin adapter + CMake detection live in this repo.

### Matrix-free Newton–Krylov (state solve)
```cpp
#include <yabi/bdf_integrator.hpp>
#include <yabi/matrix_free_solver.hpp>

auto *solver = new yabi::GmresLinearSolver<double>();
solver->setPreconditioner(yabi::KrylovPreconditioner::Jacobi); // default
solver->setTolerance(1e-12, 0.0, 200, 30);                     // rtol, atol, max it, restart
integrator.setLinearSolver(solver); // takes ownership; auto-enables matrix-free assembly
```
`FgmresLinearSolver` and `BicgstabLinearSolver` are drop-in alternatives.
The Krylov iteration drives the **state** Newton solve through the RHS
Jacobian–vector product (`AbstractRHS::jacobianVectorProduct`, one CppAD
`Forward(1)` sweep for AD tapes) and never sees a stored dense `J`. The scaled
CSC `P` is still assembled for the Newton scaling and the preconditioner
(`Jacobi` or the exact `Direct` solve). The path requires an RHS sparsity
pattern; without one the historical dense solve is used and a single warning
is logged instead of silently downgrading. `operatorApplications()` reports
how many `J·v` products ran, so a test can prove the Krylov path was taken.

Without CppAD, supply the same information on `UserODERHS` explicitly:

```cpp
#include <yabi/user_ode_rhs.hpp>

yabi::UserODERHS<double> rhs(n, f, jac);
yabi::SparsityPattern pattern;          // column-oriented df/dy structure
pattern.rows = pattern.cols = n;
pattern.column_rows = {{0, 1}, {0, 1}}; // per-column row indices
rhs.setSparsityPattern(pattern);        // validated; throws if malformed
rhs.setJacobianVectorProduct([](const double *y, const double *v, double *Jv) {
    // analytic J(y)*v; J is never formed
});
```
Both hooks are additive: without them `UserODERHS` keeps the historical dense
behaviour. This activates the same CSC sparse path used by
`CppADRHS`, in `YABI_USE_CPPAD=OFF` builds too.

### Pluggable preconditioners, J-free solves, and inexact-Newton forcing
The built-in preconditioners remain `None`, `Jacobi`, and `Direct`, but a user
object now takes precedence. Implement `yabi::Preconditioner<Scalar>` (or wrap
callbacks with `yabi::FunctionPreconditioner`) and install it on the solver:

```cpp
#include <yabi/preconditioner.hpp>

class MyAmg : public yabi::Preconditioner<double> {
  public:
    bool requiresAssembledMatrix() const override { return true; }
    bool update(const yabi::CscMatrix<double> *P,
                const yabi::MatrixFreeOperator<double> *op,
                const std::vector<double> &diagonal, unsigned int n) override {
        // Build/reuse an ILU/AMG/block preconditioner from P (or op).
        return true;
    }
    bool apply(const double *r, double *z) override {
        // z = M^{-1} r.
        return true;
    }
};

MyAmg my_amg;
solver->setUserPreconditioner(&my_amg); // borrowed; caller owns it
```

This is the adapter boundary for external packages (AMGCL, hypre, Trilinos,
block/physics-based solvers); YABI does not vendor ILU/AMG. A preconditioner
whose `requiresAssembledMatrix()` returns `false` enables a genuinely J-free
state solve: no sparsity pattern is requested and no CSC `J`/`P` is assembled.
`update()` then receives a null matrix, the live scaled `P·x` operator, and an
empty diagonal, and `D1` is estimated with one `J·v` product in direction
`2^D2` (documented as a signed row-sum approximation). A failed Krylov solve
reports failure rather than silently falling back, because no assembled matrix
exists.

An inexact-Newton forcing term can adjust the relative linear tolerance from
the nonlinear residual:

```cpp
solver->setForcingTerm(yabi::eisenstatWalkerForcing<double>());
// or any callable taking NewtonLinearizationContext<double>
```

Without a callback the fixed `setTolerance()` value is used unchanged;
`lastRelativeTolerance()` reports the value actually used.

### Batched and matrix-free sensitivity interfaces
`AbstractRHS::jacobianMatrixProduct(x, directions, columns, result)` computes
`J·V` for several column-major directions, with a correct per-column default,
one-Jacobian evaluation for `UserODERHS`, and one multi-direction CppAD
`Forward(1)` sweep. `MatrixFreeOperator::applyBlock` is the batched operator
hook, and `yabi::SensitivityLinearization` (with
`yabi::MatrixFreeSensitivityLinearization` as an adapter) names the contract a
future sparse/matrix-free sensitivity solve consumes. These are interfaces
only: the current forward/reverse/simultaneous replay still uses explicit
`J`/`H` and its numbers are unchanged.

Sensitivities are deliberately unchanged: forward/reverse/simultaneous
sensitivity still uses the explicit `J`/`H` blocks and is not matrix-free.
`tests/component/matrix_free_newton.cpp` verifies all three Krylov methods on
HIRES against the direct trajectory and checks that the sensitivity result is
unchanged; `tests/component/matrix_free_pure.cpp` does the same through
`UserODERHS` and runs in both build configurations;
`tests/component/matrix_free_jfree.cpp` exercises the no-pattern J-free path.

### Selectable Jacobian backends (T33)
How the Newton-matrix Jacobian is formed for the *plain state* solve is chosen
explicitly through `setJacobianBackend`:

```cpp
#include <yabi/jacobian_backend.hpp>

integrator.setJacobianBackend(yabi::JacobianBackend::Sparse);  // or Dense / MatrixFree / Auto
```

- `Dense` — `J`/`P` via `Forward(1)` directional columns, dense LAPACK. Always
  available; the default solver path.
- `Sparse` — `J`/`P` in CSC via `sparseJacobian()`, factorized by a CSC-capable
  solver. Requires a RHS Jacobian sparsity pattern + `supportsCsc()` solver.
- `MatrixFree` — `J·v` via `jacobianVectorProduct()`, solved by a Krylov solver
  (sparse or J-free preconditioner). Requires a `supportsMatrixFree()` solver,
  plus a pattern if its preconditioner needs an assembled matrix.
- `Auto` (default) — the historical capability-driven selection: dense, unless
  you install a UMFPACK/KLU solver (→ sparse) or a Krylov solver (→ matrix-free).

The Jacobian *source* is orthogonal and is chosen by which `AbstractRHS` you
construct: `CppADRHS` (AD) or `UserODERHS` (analytic callbacks) both feed every
representation above, so "user-analytic" is not a separate backend value.

An explicit `Dense`/`Sparse`/`MatrixFree` choice is honoured strictly: if the
RHS + installed solver cannot satisfy it, the solve **fails loudly** rather than
silently downgrading to dense. The sparse/matrix-free paths cover the plain
state Newton solve only — sensitivity solves (simultaneous seeds, or taped
forward/reverse replay) always use the explicit dense `J`/`H` path, so nothing
is silently lost. `getJacobianBackend()` reports the selection.

### Pure/manual ODE interface (no CppAD dependency in user code)
Users provide ODE as `std::function<void(double,const double*,int,double*)>` via `include/yabi/user_ode.hpp`. See `tests/system/manual_sensitivity_pure.cpp`.

## Tests
```sh
ctest --test-dir build --output-on-failure   # all tiers
ctest --test-dir build -L unit               # isolated component contracts
ctest --test-dir build -L component          # native sparse, matrix-free, SuiteSparse
ctest --test-dir build -L system             # mathematical benchmarks
./run_tests.sh                                # verbose full run
```

## Detailed usage

### 1) Providing the ODE (user function)
- **Via CppAD tape (AD):** Build `CppAD::ADFun<double>` from RHS `f(x)` (or `f(t,x)`) — see Usage example. CppAD optional.
- **Via analytical derivatives / manual interfaces:** Define RHS and provide analytic Jacobian/derivative routines through `include/yabi/user_ode.hpp`. See `examples/manual_sensitivity_pure.cpp` and `examples/manual_sensitivity_demo.cpp` (no CppAD tape required in pure form).

### 2) Analytical derivatives
- Supply analytic Jacobian `df/dy` and/or derivative RHS for sensitivity equations as appropriate for your formulation.
- Use the user-ODE interface; see `include/yabi/user_ode.hpp` for the API and `examples/manual_sensitivity_pure.cpp` for a complete example.

### 3) Pluggable linear solvers (dense + sparse)
**Interface:** `include/yabi/linear_solver.hpp` defines the abstract solver used by the corrector/linear solves. Implement `initialize()`, `factorize()` and `solve()`; a solver that consumes native CSC storage additionally overrides `supportsCsc()`, `factorizeCsc()` and `solveCsc()`.

**Built-in:**
- `DenseLUSolver<double>` (default, LAPACK `dgetrf`/`dgetrs`)
- `UmfpackLinearSolver<double>` (SuiteSparse UMFPACK)
- `KluLinearSolver<double>` (SuiteSparse KLU)
- `GmresLinearSolver<double>` / `FgmresLinearSolver<double>` / `BicgstabLinearSolver<double>`
  (matrix-free state Newton solve via `J·v`; see above)

**Connect your own solver (complete example):**
```cpp
#include <yabi/linear_solver.hpp>
#include <yabi/bdf_integrator.hpp>

class MySparseSolver : public yabi::LinearSolver<double> {
public:
  bool initialize(unsigned int n) override { /* size workspace */ return true; }

  // Dense interface: P is column-major n x n.
  bool factorize(double *P, int *ipiv, long *scalingD1, long *scalingD2,
                 unsigned int n) override { /* build symbolic/numeric */ return true; }
  int solve(double *P, int *ipiv, long *scalingD1, long *scalingD2,
            double *b, unsigned int n, unsigned int nrhs) override { /* solve */ return 0; }

  // Optional native CSC path (T31): return true to receive CscMatrix<double>.
  bool supportsCsc() const override { return true; }
  bool factorizeCsc(const yabi::CscMatrix<double>& P, int *ipiv,
                    long *scalingD1, long *scalingD2, unsigned int n) override { return true; }
  int solveCsc(const yabi::CscMatrix<double>& P, int *ipiv,
               long *scalingD1, long *scalingD2, double *b,
               unsigned int n, unsigned int nrhs) override { return 0; }
};

yabi::BDFIntegrator<double> integrator(rhs);
integrator.setLinearSolver(new MySparseSolver()); // integrator takes ownership
```
See `include/yabi/detail/dense_linear_solver.ipp`, `umfpack_linear_solver.ipp`, `klu_linear_solver.ipp` for reference implementations. Any sparse solver (MUMPS, SuperLU, Pardiso, Cholmod, etc.) can be plugged this way.

**Complete examples:** `examples/heat1d_umfpack_literature.cpp`, `examples/heat1d_dense_literature.cpp`, `examples/manual_sensitivity_pure.cpp`, `examples/manual_sensitivity_demo.cpp`.
