Metadata-Version: 2.4
Name: semiflow-pde
Version: 0.15.0b0
Classifier: License :: OSI Approved :: MIT License
Classifier: License :: OSI Approved :: Apache Software License
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Programming Language :: Python :: 3
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 :: 3.15
Classifier: Operating System :: OS Independent
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Typing :: Typed
Requires-Dist: numpy>=1.24
Summary: Python bindings for SemiFlow: evolution-equation / PDE engines via Chernoff approximation of operator semigroups
Keywords: pde,semigroup,chernoff,diffusion,numerical,simulation
Author-email: Ilia Volkov <ilia.volkov@outlook.com>
License: MIT OR Apache-2.0
Requires-Python: >=3.10, <3.16
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM
Project-URL: Changelog, https://github.com/VolkovIlia/semiflow/blob/master/CHANGELOG.md
Project-URL: Documentation, https://github.com/VolkovIlia/semiflow#readme
Project-URL: Homepage, https://github.com/VolkovIlia/semiflow
Project-URL: Issues, https://github.com/VolkovIlia/semiflow/issues
Project-URL: Release Notes, https://github.com/VolkovIlia/semiflow/releases
Project-URL: Repository, https://github.com/VolkovIlia/semiflow

<!-- Generated by `cargo xtask readme` from docs/readme/README.src.md. Do not edit by hand: edit the source and regenerate. -->
# semiflow-pde — Python bindings for SemiFlow

[![CI](https://github.com/VolkovIlia/semiflow/actions/workflows/ci.yml/badge.svg)](https://github.com/VolkovIlia/semiflow/actions/workflows/ci.yml)
[![PyPI](https://img.shields.io/pypi/v/semiflow-pde)](https://pypi.org/project/semiflow-pde/)
[![License](https://img.shields.io/badge/license-MIT%20OR%20Apache--2.0-blue)](#license)
[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.20837851.svg)](https://doi.org/10.5281/zenodo.20837851)

> **Status: beta (0.15.0-beta).** The API is stabilising toward 1.0; minor
> versions may still make breaking changes. Bug reports and feedback are welcome.

**SemiFlow solves evolution equations `∂ₜu = Lu` (heat, diffusion,
Schrödinger, graph and manifold PDEs, …) by Chernoff approximation of operator
semigroups: no matrix exponentials, no linear solves, a flat memory footprint.**
The core is a `no_std + alloc` Rust library with 3 runtime
dependencies (`num-traits`, `libm`, `num-complex`); the same engines are available from Python,
JavaScript/WebAssembly and C.

This package is the Python binding (PyO3, abi3 wheels for CPython ≥ 3.10) of
the Rust crate [`semiflow`](https://crates.io/crates/semiflow). Every class
below runs the Rust engine; the GIL is released during `evolve`.

The method evaluates `e^{tL}f` by iterating an explicit step operator `S(τ)`:
`(S(t/n))ⁿ f → e^{tL} f` as `n → ∞`. Each step is an explicit stencil over a
`Vec`-backed grid with reused scratch buffers, so steady-state evolution does
not allocate.

> Mathematical foundation: **Theorem 6 of I. D. Remizov (2025)**, *Vladikavkaz
> Math. J.* 27(4), 124–135 ([doi:10.46698/a3908-1212-5385-q](https://doi.org/10.46698/a3908-1212-5385-q)).

> **Honest performance note.** The primary measured advantage is **memory
> frugality**: a flat ~3 MB working set across the 1D/2D/3D families versus
> 50–418 MB for heavy frameworks. Wall-clock speed is **not** a general
> strength — for matched-accuracy PDE solving SemiFlow is slower than adaptive
> ODE solvers and spectral methods, and parallelism pays off only on large 3D
> grids. Two measured niches: tail-latency-sensitive pricing
> (`Diffusion4thChernoff` ≈ 41 ns p99.9 per tick) and L1-resident low-rank
> carriers (`TtChernoff`, `ReverseChernoff`). Treat memory and latency as
> measured properties of the concrete grid types, not guarantees of the traits.
> The 2026-10 overhaul ([`docs/perf/overhaul-v0_15.md`](https://github.com/VolkovIlia/semiflow/blob/master/docs/perf/overhaul-v0_15.md))
> made the grid engines 2.6–39× faster bit-for-bit, the stiff Chebyshev graph
> action `O(√(λt))` instead of `O(λt)` (559× on issue #16), ETDRK4 2.1× and the
> general Taylor action 2.1× faster; Richardson extrapolation (`richardson=`)
> raises any engine's order by `levels − 1`. Head-to-head comparisons quoted in
> the user guide predate it.

## Install

```sh
pip install semiflow-pde
```

The import name is `semiflow`. Wheels are abi3 (one wheel per platform covers
CPython 3.10+). To build from a source checkout instead (needs a Rust toolchain):

```sh
pip install maturin
maturin develop --release -m crates/semiflow-py/Cargo.toml
```

## Quickstart

Solve `∂ₜu = ∂ₓₓu` from `u₀(x) = e^{-x²}` to `t = 1` and compare with the
closed-form solution `u(1,x) = 5^{-1/2}·e^{-x²/5}`:

```python
import numpy as np
import semiflow

n = 1000
xs = np.linspace(-10.0, 10.0, n)
state = semiflow.Heat1D(-10.0, 10.0, n, np.exp(-xs**2))
state.evolve(t=1.0, n_steps=100)   # releases the GIL while Rust runs

exact = np.exp(-xs**2 / 5.0) / np.sqrt(5.0)
err = np.max(np.abs(state.values() - exact))
assert err < 5e-4, err
```

Every example on this page is executed by the project's CI against the
current build.
## Array I/O conventions

- All real-valued state arrays are `numpy.float64` (`np.float64`).
- Schrödinger and `SchrodingerComplex1D` state arrays are `numpy.complex128`.
- 2D state is flat `float64` in row-major x-fastest order: index `j*nx + i`
  corresponds to `u(x_i, y_j)`.
- 3D state is flat x-fastest: index `k*nx*ny + j*nx + i`.
- `values()` always returns a **copy** of the internal Rust state; mutations
  to the returned array do not affect the object.
- Inputs are validated for `NaN`/`Inf` at construction and before `evolve`;
  non-finite inputs raise `SemiflowError(kind='NanInf')`.
- All finite-check and grid-size errors raise `SemiflowError`.

## Error model

All semiflow-py operations raise a single exception type:

```python
from semiflow import SemiflowError
```

The `.kind` attribute (a string) identifies the error category:

| `kind` | When raised |
|--------|-------------|
| `GridMismatch` | Invalid geometry, mismatched array lengths |
| `NanInf` | Input array contains NaN or Inf |
| `OutOfDomain` | Parameter out of valid range (e.g. `t < 0`, `n < 4`), or an unrecognised value for most string keywords (e.g. `boundary=`, `dtype=`, `kernel=`, `method=`, `phi=`, `solver=`, `precond=`, `reaction_solver=`, `nonlinearity=`) |
| `BoundaryFailure` | Reserved: grid under-resolved for the Chernoff shift (core `GridUnderresolved`); not raised by the current API |
| `CflViolated` | CFL-like stability constraint exceeded |
| `ConvergenceFailed` | Magnus / adaptive integration check failed, or an iteration cap was hit (PCG in `SpdResolvent` and `path="implicit"`, Arnoldi substeps of `method="krylov"` / `phi="krylov"`, `RationalKrylov` / `path="rational"` needing more than `m_max`, Rosenbrock steps of `reaction_solver="rosenbrock"`) |
| `Unsupported` | Unrecognised `path=`, `subordinator=`, `manifold=` or reaction-model name, or a combination that is not implemented (e.g. `solver="tridiagonal"` on a non-tridiagonal operator, `path="implicit"` on `GraphKrylov`) |
| `Panic` | Unrecoverable internal Rust panic (should never occur) |

## Boundary policies

All 1D/2D/3D kernels accept a keyword argument `boundary`:

| Value | Semantics | Typical use |
|-------|-----------|-------------|
| `"reflect"` (default) | Mirror / zero-flux Neumann at grid boundaries | General PDEs; required by G1/G2 oracle tests |
| `"periodic"` | Periodic wrap with period `n·dx` (node `n` is node `0`: sample `[a, b)` with `b = a + n·dx`) | Periodic domains |
| `"zero"` | Extend with 0.0 outside domain | Solutions that vanish at the boundary (barriers, puts far OTM in log-space) |
| `"linear"` | Linear extrapolation from the two outermost nodes | **Asymptotically-linear far-field** (European calls, linear ramps) |

### Far-field / Dirichlet-like boundaries for finance

Users coming from classical finite-difference option pricers often look for an
inhomogeneous Dirichlet far-field ("set `V = S − K e^{−rτ}` at `S_max`") and,
finding only `"zero"`, conclude the library cannot price calls.  In fact
`boundary="linear"` is the correct and recommended closure for European call
pricing.

**Why `"linear"` works for calls:**  The Black-Scholes / CEV solution satisfies
`V(S, τ) ≈ S − K e^{−rτ}` as `S → ∞`.  This asymptote is linear in `S`, so
`V_SS → 0`.  Linear extrapolation from the two outermost grid nodes reproduces
an affine far-field exactly (to machine precision), introducing no spurious
curvature or kink.  Validated on `Shift1D.with_arrays` with `S ∈ [0, 4K]`,
`n = 1025`: ATM relative error ≈ **8.5e-5** with no boundary artifact.

**Puts in log-price space:** When a European put is priced in log-price
coordinates (`z = ln S`), the solution satisfies `u_z → 0` at both ends of a
sufficiently wide domain.  Either `"reflect"` (zero-flux Neumann) or `"linear"`
works; `"reflect"` is marginally preferable because it imposes the correct
zero-derivative condition exactly.

**Example — Black-Scholes call via `Shift1D`:**

```python
import numpy as np
import semiflow as rp

# Black-Scholes PDE: dV/dt = (sigma^2 / 2) S^2 V_SS + r S V_S - r V
# Rewrite as L V = a(S) V_SS + b(S) V_S + c(S) V with:
#   a(S) = 0.5 * sigma^2 * S^2,  b(S) = r * S,  c(S) = -r

K, T, r, sigma = 100.0, 1.0, 0.05, 0.20
S_max = 4.0 * K      # wide enough that far-field is linear
n = 1025
S_nodes = np.linspace(0.0, S_max, n)

a_arr = 0.5 * sigma**2 * S_nodes**2
b_arr = r * S_nodes
c_arr = np.full(n, -r)

# European call payoff at maturity
u0 = np.maximum(S_nodes - K, 0.0)

state = rp.Shift1D.with_arrays(
    0.0, S_max, n,
    a_arr, b_arr, c_arr,
    c_norm_bound=r,  # upper bound on |c(x)|
    u0=u0,
    boundary="linear",   # <-- correct far-field closure for calls
)
state.evolve(t=T, n_steps=200)

V = state.values()
S_atm_idx = np.argmin(np.abs(S_nodes - K))
print(f"ATM call price (semiflow): {V[S_atm_idx]:.6f}")
# rel. error vs Black-Scholes closed form ≈ 8.5e-5 at n=1025, n_steps=200
```

**Inhomogeneous Dirichlet (future work):** The Rust core provides
`BoundaryPolicy::Dirichlet { value }` (a constant ghost-node extension), but
it is not yet wired into the Python `boundary=` string parser.  If you need to
pin `u = g_lo` at the left edge and `u = g_hi` at the right edge with distinct
values, use `"linear"` as the near-exact idiom for asymptotically-linear
payoffs.  Support for `boundary=("dirichlet", value)` is tracked in issue #20.

---

## Usage examples

### 1. `SchrodingerComplex1D` — native complex128 wavefunction

Solve `i ψ_t = (−½∂²_x + V) ψ` and verify unitarity:

```python
import numpy as np
import semiflow as rp

n = 512
xs = np.linspace(-10.0, 10.0, n)
psi0 = np.exp(-xs**2 / 2.0).astype(np.complex128)  # normalised Gaussian
dx = xs[1] - xs[0]
psi0 /= np.sqrt(np.sum(np.abs(psi0)**2) * dx)       # L2-normalise

sch = rp.SchrodingerComplex1D(-10.0, 10.0, n, psi0)
norm0 = sch.norm_squared()

sch.evolve(t=0.5, n_steps=200)

psi_t = sch.values()    # complex128 ndarray
assert abs(sch.norm_squared() / norm0 - 1.0) < 1e-12, "unitarity violated"
print(f"norm ratio = {sch.norm_squared() / norm0:.15f}")
```

### 2. `Manifold2D` — Riemannian manifold heat kernel

Solve `∂_t u = Δ_{S²} u` on the 2-sphere via MMRS 2023 Chernoff formula:

```python
import numpy as np
import semiflow as rp

nx, ny = 32, 64
u0 = np.zeros(nx * ny, dtype=np.float64)
u0[nx * (ny // 2) + nx // 2] = 1.0  # delta-like at chart centre

sphere = rp.Manifold2D(
    0.1, np.pi - 0.1, nx,    # theta axis
    0.0, 2 * np.pi,   ny,    # phi axis
    u0,
    manifold="sphere2",
    radius=1.0,
    curvature_correction=True,  # enables R/12 correction -> order 2
)
sphere.evolve(t=0.02, n_steps=50)

u_t = sphere.values()    # float64 ndarray, length nx*ny (row-major theta-fastest)
print(f"integral ≈ {u_t.sum() * (np.pi / nx) * (2 * np.pi / ny):.4f}")
```

Available manifolds: `"torus"` (flat T²), `"sphere2"` (S²(r)), `"hyperbolic2"`
(Poincaré disk H²(s)). The `radius` parameter sets r or s.

### 3. `ReactionDiffusion2D` — Gray–Scott patterns, or your own kinetics

Any number of species, each with its own diffusivity; built-in kinetics by
name (`fisher_kpp`, `allen_cahn`, `nagumo`, `gray_scott`, `fitzhugh_nagumo`,
`brusselator`, `linear`) or any NumPy-vectorised callable `f(t, x, u) -> du`,
called once per Runge–Kutta stage with all nodes (`x`: `(dim, N)`, `u`:
`(K, N)`), never once per node:

```python
import numpy as np
import semiflow

n = 128
# periodic: node n is node 0, so the unit torus is sampled on [0, 1)
x, y = np.meshgrid(np.linspace(0, 1, n, endpoint=False), np.linspace(0, 1, n, endpoint=False))
bump = np.exp(-((x - 0.5)**2 + (y - 0.5)**2) / 0.005)
u0 = np.stack([1.0 - 0.5 * bump, 0.25 * bump])          # shape (K=2, ny, nx)

rd = semiflow.ReactionDiffusion2D(
    0.0, 1.0 - 1.0 / n, n, 0.0, 1.0 - 1.0 / n, n, u0,
    diffusivity=[2e-5, 1e-5], reaction="gray_scott",
    params={"feed": 0.037, "kill": 0.06}, boundary="periodic",
)
rd.evolve(2000.0, n_steps=2000)            # Strang splitting, order 2
u, v = rd.values()

# Custom kinetics + a source term, with Richardson extrapolation (order 3):
def kinetics(t, xs, u):
    return np.stack([u[0] * (1 - u[0]) - u[0] * u[1] + 0.1 * np.sin(t),
                     0.5 * u[0] * u[1] - 0.2 * u[1]])

rd2 = semiflow.ReactionDiffusion1D(-10, 10, 401, np.full((2, 401), 0.5),
                                   diffusivity=[1.0, 0.1], reaction=kinetics)
rd2.evolve(5.0, n_steps=50, richardson=2)

# Stiff chemistry (Robertson): adaptive L-stable Rosenbrock reaction steps.
def robertson(t, xs, u):
    a, b, c = u
    return np.stack([-0.04 * a + 1e4 * b * c,
                     0.04 * a - 1e4 * b * c - 3e7 * b * b, 3e7 * b * b])

u0 = np.stack([np.ones(8), np.zeros(8), np.zeros(8)])
rd3 = semiflow.ReactionDiffusion1D(0, 1, 8, u0, diffusivity=[0.0, 0.0, 0.0],
                                   reaction=robertson, reaction_solver="rosenbrock",
                                   rtol=1e-5, atol=1e-12)
rd3.evolve(40.0, n_steps=1)   # RK4 needs ~25 000 substeps per half-step here (NaN at <= 21 000)
assert abs(rd3.values()[0, 0] - 0.7158270687) < 1e-4

# Cross-diffusion: a K×K diffusivity matrix D (diagonalisable, real eigenvalues
# >= 0, verified); each diffusion step diffuses the modes of D exactly.
rd4 = semiflow.ReactionDiffusion1D(0, 1, 64, np.full((2, 64), 0.5),
                                   diffusivity=[[1e-2, 5e-3], [0.0, 2e-3]],
                                   reaction="gray_scott", params={"feed": 0.04, "kill": 0.06})
rd4.evolve(1.0, n_steps=20)
```

---

## Class reference

Classes are grouped by kernel family. All stateful classes expose at least
`evolve(t, n_steps=100)` (mutates in-place, GIL released) and `values()` →
`NDArray[np.float64]` (copy). See `__init__.pyi` for complete signatures.

### 1D diffusion family

| Class | Kernel | Order | Notes |
|-------|--------|-------|-------|
| `Heat1D` | `DiffusionChernoff` | 2 (`1 + richardson`) | Unit or variable-`a`; `.with_a_array` / `.with_a_function` factories; `.evolve(t, n_steps, richardson=L)` Richardson extrapolation (`1 ≤ L ≤ 6`, f64 only, ADR-0207) |
| `Heat1D4th` | `Diffusion4thChernoff` | 4 | 4th-order temporal; `.with_a_array` |
| `Heat1D6th` | `Diffusion6thChernoff` | 6 | 6th-order temporal; `.with_a_array` |
| `Heat1DZeta4` | `Diffusion4thZeta4Chernoff` | 4 | ζ⁴ kernel; `.with_quintic_sampling()` opt-in |
| `Heat1DZeta6` | `Diffusion6thZeta6Chernoff` | 6 | ζ⁶ kernel; Quintic spatial unconditional |
| `Heat1DZeta8` | `Diffusion8thZeta8Chernoff` | 8 | ζ⁸ kernel; Chebyshev sampling default |
| `TruncatedExp1D` | `TruncatedExpDiffusionChernoff` | 2 | CFL-conditional truncated-exp |
| `TruncatedExp4th1D` | `TruncatedExp4thDiffusionChernoff` | 4 | 4th-order truncated-exp |
| `DriftReaction1D` | `DriftReactionChernoff` | 2 | `b(x) ∂_x u + c(x) u`; `.with_arrays` |
| `Shift1D` | `ShiftChernoff1D` | 1 | Universal `a ∂² + b ∂ + c`; `.with_arrays` |
| `Strang1D` | `StrangSplit` (diffusion + drift) | 2 | Advection-diffusion `∂²u + b ∂u`; default `b=0.5` |

### Reaction–diffusion systems (ADR-0208)

| Class | Kernel | Order | Notes |
|-------|--------|-------|-------|
| `ReactionDiffusion1D` | `ReactionDiffusion` + `DiffusionChernoff` | 2 (`1 + richardson`) | `K` species, `u0` `(K, n)`; built-in kinetics or callable `f(t, x, u)`; `reaction_solver="rosenbrock"` (`rtol`, `atol`, `jacobian=`) for stiff kinetics (ADR-0212); `K×K` (`K ≤ 16`) `diffusivity` matrix for cross-diffusion (ADR-0213) |
| `ReactionDiffusion2D` | … + `Strang2D` | 2 | `u0` `(K, ny, nx)` |
| `ReactionDiffusion3D` | … + `Strang3D` | 2 | `u0` `(K, nz, ny, nx)` |

### Operator splitting — multi-dimensional

| Class | Kernel | Order | Notes |
|-------|--------|-------|-------|
| `Heat2D` | `Strang2D` | 2 | Unit diffusion on 2D grid; flat x-fastest output |
| `Heat3D` | `Strang3D` | 2 | Unit diffusion on 3D grid; flat x-fastest output |
| `Heat2DVarA` | `Strang2D` + variable-a | 1 | **Non-divergence** `a_x(x) u_xx + a_y(y) u_yy`; pass `a_x`, `a_y` arrays. Order 1, not 2: the axis kernels freeze `a` at the node (ADR-0191 AM1) |
| `Heat3DVarA` | `Strang3D` + variable-a | 1 | **Non-divergence** `a_x u_xx + a_y u_yy + a_z u_zz`; same order caveat as `Heat2DVarA` |
| `NonSeparable2D` | 5-leg palindromic | 2 | `∂²_x + ∂²_y + c·∂_x ∂_y`; scalar or `.with_beta_array` |
| `NonSeparable2DAniso` | 5-leg + position-dep. β | 2 | `∂²_x + ∂²_y + β(x,y)·∂_x ∂_y`; requires `beta_values` array |

### Schrödinger

| Class | Kernel | Notes |
|-------|--------|-------|
| `Schrodinger1D` | `SchrodingerChernoff<f64>` | Real-pair split; `values()` → `complex128` |
| `SchrodingerComplex1D` | `SchrödingerChernoffComplex` | Native `complex128` state; exact unitary (ADR-0079 Option B) |

Both support `.with_potential(v_array)` and `.norm_squared()`.

### Boundary-condition kernels

| Class | Kernel | Order | Physics |
|-------|--------|-------|---------|
| `Resolvent1D` | `LaplaceChernoffResolvent` | — | `(λI − ∂²)⁻¹ g` via GL-32 quadrature; `.eval(lambda_, g)` + `.residual(lambda_, g)` |
| `Killing1D` | `KillingChernoff` | 1 | Absorbing (Dirichlet) BC via Feynman-Kac; `lo`/`hi` kwargs |
| `Reflected1D` | `ReflectedHeatChernoff` | 2 | Neumann (reflecting) BC via Walsh 1986 image method; `origin` kwarg |
| `Robin1D` | `RobinHeatChernoff` | 1 | Robin BC `α u − β ∂_n u = 0`; `alpha`, `beta`, `origin` kwargs |

### Time-dependent and subordinated

| Class | Kernel | Notes |
|-------|--------|-------|
| `Howland1D` | `HowlandLift<DiffusionChernoff>` | Nonautonomous lift (Howland 1974); `n_t`, `t_horizon` kwargs; `.evolve()` takes no args |
| `Subordinated1D` | `SubordinatedChernoff` | Bochner-Phillips subordination (Butko 2018); backends: `"stable"`, `"gamma"`, `"inverse_gaussian"` |

### Geometry and hypoelliptic operators

| Class | Manifold / Group | Notes |
|-------|-----------------|-------|
| `Manifold2D` | Torus / S²(r) / H²(s) | MMRS 2023 formula with optional R/12 correction; `manifold=`, `radius=`, `curvature_correction=` kwargs |
| `HypoellipticChernoffKolmogorov` | Kolmogorov phase space | `∂_t p = v ∂_x p + ½ ∂²_v p`; 2D state `nx×nv` |
| `HypoellipticChernoffEngel` | Engel step-3 Carnot (ℝ⁴) | `n**4` flat state; `n` per-axis |
| `HypoellipticChernoffHeisenberg` | Heisenberg H₁ | `.kernel(h, x, y, tc)` point evaluator; `heisenberg_heat_kernel(h, x, y, tc)` free function |

### Graph PDE

| Class / Function | Role |
|-----------------|------|
| `Graph.path(n)` / `.cycle(n)` / `.from_edges(n, edges)` / `.erdos_renyi(n, p, seed)` | Graph topology builders |
| `GraphPath(n)` | Legacy path builder (use `Graph.path(n)`) |
| `Laplacian.combinatorial(graph)` / `.normalized(graph)` | Laplacian assembly |
| `GraphHeat(graph=..., laplacian=..., rho_bar=...)` | Order-1 static graph heat |
| `GraphHeat4th(graph=..., laplacian=..., rho_bar=...)` | Order-4 static |
| `GraphHeat6(graph=..., laplacian=..., rho_bar=...)` | Order-6 static |
| `MagnusGraphHeat(graph, lap_at_t, rho_bar)` | Magnus K=4 time-varying |
| `MagnusGraphHeat6(graph=..., laplacian=..., lap_at_t=..., rho_bar_max=...)` | Magnus K=6 |
| `VarCoefGraphHeat(graph, a, rho_bar)` | Variable node-conductivity |
| `VarCoefMagnusGraph(n_nodes, lap_at_t=..., a_at_t=..., rho_bar_max=..., a_sup_max=...)` | Variable-coef Magnus K=4 |
| `QuantumGraph.path(n_edges)` / `.star(n_arms)` / `.from_edges(edges)` | Metric graph (edge lengths) |
| `QuantumGraphHeat(qgraph)` | Kirchhoff-vertex heat Chernoff |
| `GraphTraj(graph, t_horizon)` | Fixed-topology graph trajectory |
| `StrangGraph.from_path(graph)` / `.from_cycle(graph)` | Palindromic Strang split on graph |

### Matrix and point-eval kernels

| Class / Function | Role |
|-----------------|------|
| `MatrixDiffusion1D(xmin, xmax, n, u0, *, a_diag, c_coupling)` | Coupled 2-component 1D diffusion; flat state length `2*n` |
| `PointEval(xmin, xmax, n)` | Pointwise evaluation via Backend A; `.eval_at(tau, u0, x, n_steps)` |
| `sample_gridfn2d(values, x0min, x0max, nx, x1min, x1max, ny, cx, cy)` | Bilinear interpolation at chart point |

### Anisotropic multi-D

| Class | Notes |
|-------|-------|
| `AnisotropicShiftND2(nx, ny, xmin, xmax, ymin, ymax, a_values, *, b_values, c_values)` | 2D anisotropic shift; order 1 (ADR-0112); `a_values` is flat `2×2×nx×ny` SPD tensor |
| `AnisotropicShiftND3(nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax, a_values, *, b_values, c_values)` | 3D variant |

### Adjoint and adaptive wrappers

| Class | Notes |
|-------|-------|
| `Adjoint(xmin, xmax, n, u0, *, kernel="heat2", self_adjoint=False, boundary="reflect")` | Adjoint semigroup; `kernel` in `"heat2"`, `"heat4"`, `"heat6"`, `"drift"`, `"shift"` |
| `AdaptivePI(xmin, xmax, n, u0, *, kernel="heat2", tol_abs=1e-6, tol_rel=1e-4, boundary="reflect")` | PI-controller adaptive step |

### Reverse-mode AD (ADR-0156)

| Class | Notes |
|-------|-------|
| `ReverseHeat1D(theta, xmin, xmax, n_grid, n_steps)` | Reverse-mode AD for constant-a 1D heat (narrow scope: constant-a `DiffusionChernoff` only, §51.5); `.value_and_grad(tau, u0, target) -> (float, float)` |

**Constructor parameters:**

| Parameter | Type | Constraint |
|-----------|------|------------|
| `theta` | `float` | Diffusivity θ > 0, finite |
| `xmin` | `float` | Left domain boundary |
| `xmax` | `float` | Right domain boundary (xmax > xmin) |
| `n_grid` | `int` | Grid nodes (>= 4) |
| `n_steps` | `int` | Chernoff steps per `.value_and_grad` call (>= 1) |

**`.value_and_grad(tau, u0, target) -> (float, float)`:**

| Parameter | Type | Notes |
|-----------|------|-------|
| `tau` | `float` | Per-step time increment (> 0, finite) |
| `u0` | `numpy.ndarray[float64]` | Initial condition, length `n_grid` |
| `target` | `numpy.ndarray[float64]` | Target state, length `n_grid` |
| returns `value` | `float` | L² loss `‖(F_θ(τ))ⁿ u₀ − target‖²` |
| returns `grad` | `float` | `∂J/∂θ` (K=1 forward-mode Dual; 0-ULP vs core, §51.4) |

```python
import numpy as np
import semiflow as rp

n_grid = 24
xs = np.linspace(-4.0, 4.0, n_grid)

rc = rp.ReverseHeat1D(theta=0.4, xmin=-4.0, xmax=4.0, n_grid=n_grid, n_steps=8)
u0     = np.exp(-xs**2)
target = np.zeros(n_grid)

value, grad = rc.value_and_grad(tau=0.05, u0=u0, target=target)
print(f"loss={value:.6e}  ∂J/∂θ={grad:.6e}")
```

Raises `SemiflowError` with `.kind` in `{'OutOfDomain', 'GridMismatch', 'NanInf'}`.

**NARROW scope (§51.5):** constant-a `DiffusionChernoff` only; θ is the
uniform diffusivity. Variable-coefficient and nonlinear kernels are out of scope
for `ReverseHeat1D`.

### Tensor-train carriers (ADR-0171 / ADR-0178)

Curse-escaped O(d·n·r²) storage for separable diagonal-A diffusion on ℝᵈ.
State and evolvers share the `TtState` carrier; `TtEvolver` / `TtCoupledEvolver`
cover constant-coefficient cases, `VarCoefTtEvolver` covers variable-coefficient
additive-separable operators (issue #2, ADR-0178).

| Class | Constructor | Key method | Notes |
|-------|-------------|------------|-------|
| `TtState(slices)` | `slices: list[NDArray[np.float64]]` — per-axis 1-D arrays | `.ndim()`, `.n_j(j)`, `.peak_rank()`, `.storage_size()`, `.inner_separable(functionals)` | Shared carrier for all TT evolvers |
| `TtEvolver(a, b, c, dom_min, dom_max, eps_round)` | `a`, `b`: `list[float]` per-axis coeffs; `c: float`; bounds lists; `eps_round: float` | `.evolve(state, t_final, n_steps)` — mutates `TtState` in-place | Diagonal-A Gaussian class (§52) |
| `VarCoefTtEvolver(a_axis, b_axis, v_axis, domain, eps_round)` | `a_axis`, `b_axis`, `v_axis`: `list[list[float]]` per-axis nodal values; `domain: list[tuple[float,float]]`; `eps_round: float` | `.evolve(state, t_final, n_steps)` | Variable-coef additive-separable; rank-1 IC → rank-1 output (ADR-0178) |
| `TtCoupledEvolver(a, b, c, coupling, dom_min, dom_max, eps_round)` | Same as `TtEvolver` + `coupling: tuple` — `("None",)`, `("Tridiagonal", rho)`, or `("Pairs", [(j,k,rho),...])` | `.evolve(state, t_final, n_steps)` | Nearest-neighbour pair-factor coupling (§52.9) |

### Gridless / particle carriers (ADR-0171)

Sparse signed-measure ensemble on ℝ (D=1). Curse-escaped: the 3ᴰ dense tree is
never materialised; only sparse marginals and scalar observables cross the
Python boundary.

| Class | Constructor | Key method | Notes |
|-------|-------------|------------|-------|
| `MeasureState(positions, weights, dim)` | `positions`, `weights`: `NDArray[np.float64]`; `dim: int` (must be 1) | `.n_diracs()`, `.total_variation()`, `.second_moment()`, `.marginal(axis) -> (pos, wgt)` | Particle carrier; signed-weight Dirac ensemble (§50) |
| `GridlessEvolver(a, b, c, *, voronoi_cap=64, gaussian_background=False)` | `a`, `b`, `c`: `float`; optional particle-cap / background kwargs | `.evolve(state, t_final, n_steps)` — mutates `MeasureState`; `.apply(tau, src, dst)` one-step | 1-D 3-branch Chernoff kernel; `WeightedVoronoi` reduction (ADR-0155) |

### v3 Evolver surface

| Class | Notes |
|-------|-------|
| `EvolverHeat1DUnitV3(domain_lo, domain_hi, n_grid, u0, n_chernoff)` | Zero-alloc `apply_into` hot path; `.evolve_into(t, buf)` |
| `GrowthV3` | Growth bound `(multiplier, omega)` returned by `.growth()` |

### Greeks and hyper-dual AD

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `EvolverHeat1DGreeksV3` | `(domain_lo, domain_hi, n_grid, u0, n_chernoff, scale_theta=0.5)` | `.greeks(t) -> (value, delta, gamma)`, `.size()`, `.n_chernoff()` | Hyper-dual `Dual<Dual<f64>>` sweep; returns primal + ∂u/∂θ + ∂²u/∂θ² in one pass (ADR-0133 A3) |

### Adjoint Fokker-Planck (particle pushforward)

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `AdjointFokkerPlanckV8` | `(a, b, c)` — constant diffusion, drift, reaction coefficients | `.step(tau, positions, weights, n_steps=1) -> (positions, weights)`, `.total_variation(...)`, `.second_moment(...)` | Weak-* Fokker-Planck pushforward: each Dirac δ_x spawns 4 children per step (Lemma A.1, §38.3); D=1 constant-coefficient only |

### Killing and absorbing boundary conditions

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `Killing2nd1D` | `(xmin, xmax, n, u0, *, kappa=0.0, boundary="reflect")` | `.evolve(t, n_steps=100)`, `.values()`, `.order()`, `len()` | Order-2 soft-killing `∂_t u = ∂²u − κu`; palindromic Strang split `e^{−τκ/2} C(τ) e^{−τκ/2}` (ADR-0126, §21.8) |
| `KilledDirichlet1D` | `(domain_lo, domain_hi, n_grid, u0, n_chernoff)` | `.apply(t) -> NDArray`, `.size()` | Crank-Nicolson Cayley map; absorbing Dirichlet u\|∂R=0; order 2 (ADR-0135 Amdt 2, §44.ter) |
| `DirichletHeat2nd1D` | `(xmin, xmax, n, u0, *, origin=nan, boundary="reflect")` | `.evolve(t, n_steps=100)`, `.values()`, `.order()`, `len()` | Order-2 Dirichlet via odd-image method; sibling of `Reflected1D` (Neumann); higher-order than `Killing1D` (ADR-0176, §21.9) |

### Obstacle / variational-inequality family

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `ObstacleChernoff` | `(xmin, xmax, n, u0, *, a=1.0, b=0.0, c=0.0, level=nan, obstacle_array=None)` | `.evolve(t, n_steps=100) -> NDArray`, `.values()`, `.evolve_active_set_adjoint(w_fwd, lam, tau)`, `.order()`, `len()` | `V^{n+1} = Π_g(S(Δτ)Vⁿ)`; constant or array obstacle floor; order 1 globally (§44, Theorem 44.1) |
| `ObstacleGammaV8` | `(domain_lo, domain_hi, n_grid, *, level=..., obstacle_array=None)` | `.inactive_gamma(v) -> (gamma, defined, count)`, `.size()` | Inactive-set Γ = V″ primitive; `defined[i]=False` means Γ **refused** (active set / contact); D=1 only (ADR-0153 §4.1) |
| `ObstacleNDV8` | `(xmin, xmax, nx, ymin, ymax, ny, level)` | `.apply(tau, v) -> NDArray`, `.shape()` | D=2 projective-splitting obstacle `Π_g ∘ S(Δτ)`; input shape `(nx, ny)` accepted (raveled order='F' internally); output is a flat nx*ny float64 array — use `out.reshape((nx, ny), order='F')` to recover 2D layout (ADR-0153 §4.2) |

### Resolvent time-jump family (TWS contour)

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `ResolventJumpV8` | `(domain_lo, domain_hi, n_grid, m_nodes)` | `.jump(t, g) -> NDArray`, `.size()`, `.m_nodes()` | Evaluates `e^{tA}g` via TWS parabolic-contour ILT; large-t alternative to many Chernoff steps; `m_nodes >= 6` (ADR-0138, §47) |
| `ResolventJump2DV8` | `(xmin, xmax, nx, ymin, ymax, ny, m_nodes)` | `.jump(t, g) -> NDArray`, `.shape()`, `.m_nodes()` | 2D TWS contour ILT; input/output shape `(nx, ny)` Fortran-order (ADR-0153, §47.8) |
| `ResolventJump3DV8` | `(xmin, xmax, nx, ymin, ymax, ny, zmin, zmax, nz, m_nodes)` | `.jump(t, g) -> NDArray`, `.shape()`, `.m_nodes()` | 3D TWS contour ILT; input/output shape `(nx, ny, nz)` Fortran-order (ADR-0153, §47.8) |

### Matrix diffusion (coupled 2-component)

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `MatrixDiffusion2D` | `(xmin, xmax, nx, ymin, ymax, ny, u0, *, a_diag=1.0, c_coupling=0.0)` | `.evolve(t, n_steps=100)`, `.values()`, `.order()` | 2-component coupled 2D diffusion `∂_t u = aΔu + Cu`, `C = [[0, c], [c, 0]]` (`c = c_coupling`); flat buffer length `2·nx·ny`; index `(j·nx+i)·2+component` (ADR-0124, §33.2) |
| `MatrixDiffusion3D` | `(xmin, xmax, nx, ymin, ymax, ny, zmin, zmax, nz, u0, *, a_diag=1.0, c_coupling=0.0)` | `.evolve(t, n_steps=100)`, `.values()`, `.order()` | 2-component coupled 3D diffusion `∂_t u = aΔu + Cu`, `C = [[0, c], [c, 0]]`; flat buffer length `2·nx·ny·nz`; index `(k·nx·ny+j·nx+i)·2+component` (ADR-0124, §33.3) |

### Wentzell dynamic boundary conditions

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `GammaFamily` | Static factories: `.constant(c)`, `.linear(a, b)`, `.exponential(rate)` | — | Ergonomic γ-schedule builder for `WentzellV8`; expands to pre-sampled float64 array (ADR-0153) |
| `WentzellV8` | `(domain_lo, domain_hi, n_grid, u0, n_steps, c_reaction, gamma_schedule)` or `.from_family(...)` | `.evolve(t, t_offset=0.0) -> NDArray`, `.size()`, `.n_steps()` | Dynamic Wentzell BC `∂_t u + γ(t)·∂_ν u + c·u = 0`; bulk-boundary Cayley Lie split; order 1; 1D half-line only (ADR-0151, §49) |

### Additional 1D kernels

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `DiffusionExpmv1D` | `(xmin, xmax, n, u0, *, boundary="reflect")` | `.evolve(t, n_steps=100)`, `.values()`, `.order()`, `len()` | Tolerance-driven expmv kernel (Al-Mohy & Higham 2011); `order()` returns `u32::MAX`; not a fixed convergence order (ADR-0121) |
| `DriftReaction4th1D` | `(xmin, xmax, n, u0, *, boundary="reflect")` | `.evolve(t, n_steps=100)`, `.values()`, `.order()`, `len()` | Order-4 `b(x)∂_x u + c(x)u` via palindromic `R_sym(τ/2) ∘ K5(τ) ∘ R_sym(τ/2)`; defaults `b=0.5, b'=0, c=0` (ADR-0127) |

### Sparse-grid and high-dimensional kernels

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `SmolyakD6V8` | `(domain_lo, domain_hi, n_per_axis)` | `.apply(tau, u0, n_steps=1) -> NDArray`, `.n_nodes()`, `.level()`, `.size()` | D=6 Smolyak sparse-grid unit-diffusion; level ℓ=D+3=9 (533 nodes); input/output flat `n_per_axis^6` (ADR-0138, ADR-0123 Amdt 1) |
| `ComplexTripleJumpV8` | `(domain_lo, domain_hi, n_per_axis)` | `.apply_real(tau, u0) -> NDArray`, `.verify_gamma_star() -> bool`, `.size()` | Order-4 complex triple-jump on filiform-N5 Carnot (D=5); returns real projection `Re(Ψ(τ)f)`; complex substeps are internal (ADR-0138, ADR-0136 Amdt 2) |

### Graph adjoint and Krylov families

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `GraphAdjoint` | `(graph=None, laplacian=None, *, lap_at_t, rho_bar, a=None, kernel="magnus_graph", convergence_check=True)` | `.evolve_state_adjoint(lambda_n, t, n_steps=100) -> NDArray`, `.n_nodes()` | Backward costate sweep `λ_0 = S⋆_1⋯S⋆_n · λ_n` via transpose Magnus K=4 map; supports `"magnus_graph"` / `"varcoef_magnus_graph"` (ADR-0115, §42) |
| `GraphAdjointPresampled` | `.from_presampled(graph, lap_at_t, rho_bar, n_steps, t_horizon, a=None, kernel="magnus_graph", convergence_check=True)` | `.evolve_state_adjoint(lambda_n, n_steps=None) -> NDArray`, `.evolve_state_adjoint_batched(lambda_cols, n_steps=None) -> NDArray`, `.n_nodes()`, `.n_steps()` | Pre-samples callbacks once at construction; adjoint sweep runs fully in `py.detach` with no per-step GIL re-entry (ADR-0180) |
| `GraphKrylov` | `(laplacian, *, path="chebyshev", tol=1e-10, m_max=18)` | `.evolve_batched(t, features_nc) -> NDArray`, `.n_nodes()` | Depth-independent `e^{-tL_G}·V` via Krylov; accepts `[N, C]` feature matrix; single GIL release (ADR-0185, §54) |

### Symmetric operator and FEM utilities

| Class | Constructor | Key methods | Notes |
|-------|-------------|-------------|-------|
| `SymmetricOperator` | `.from_csr(indptr, indices, data, n, sym_tol=1e-10)` | `.evolve_batched(t, v_nc, path="chebyshev", tol=1e-10, m_max=None, n_steps=100, cg_max_iter=None) -> NDArray` (`path`: `"chebyshev"`, `"lanczos"`, `"implicit"`, `"rational"`), `.n()`, `.lambda_max_bound()`, `.with_diagonal(c)`, `.to_csr()`, `.lumped_congruence(masses)`, `.resolvent(lam=0.0, mass=None, solver="auto", precond="ic0", tol=1e-12, max_iter=None) -> SpdResolvent`, `.rational_krylov(t, mass=None, tol=1e-10, gamma=None, m_max=64, solver="auto", precond="ic0") -> RationalKrylov` | Externally-assembled symmetric PSD sparse operator from CSR arrays; feeds Krylov expmv and Fréchet VJP (`symmetric_op_expmv_frechet`) (§55); `with_diagonal` adds a local reaction `A + diag(c)`, `to_csr` round-trips with `from_csr`, `resolvent` factors `(λM + A)` for steady solves (§62) |
| `RationalKrylov` | `SymmetricOperator.rational_krylov(...)` | `.apply(v) -> (NDArray, dict)`, `.apply_batched(v_nc) -> NDArray`, `.t`, `.gamma`, `.n` | `e^{-t M⁻¹A}v` by shift-and-invert Lanczos on one factor of `M/γ + A`: the solve count is an a-priori bound independent of `λ_max·t` (27 solves at `tol = 1e-10`); tridiagonal M-matrices use a cancellation-free factor (no rounding floor up to `λ_max·t = 1e10` on the gate carriers), other patterns PCG (ADR-0206, §66) |
| `ConservativeDiffusionChernoff` | `.from_k_array(n, x_lo, x_hi, k_nodes, r_contact=None, boundary="neumann")` | `.to_symmetric_operator() -> SymmetricOperator`, `.n()`, `.dx()` | Order-2 FV divergence-form `∂_x(k(x)∂_x u)` with harmonic-mean face conductivities; bridge to `SymmetricOperator` Krylov path (§56) |
| `GeneralOperator` | `.from_csr(n, indptr, indices, data)` | `.evolve_batched(t, v_nc, method="taylor", tol=1e-10, m=30) -> NDArray`, `.apply_transpose(v) -> NDArray`, `.n()`, `.norm_inf_bound()` | Externally-assembled **possibly non-symmetric** CSR operator; `e^{-tA}v` via the symmetry-agnostic Al-Mohy–Higham Taylor `expmv` (cost `Θ(t‖A‖_∞)` matvecs, §57, ADR-0195) or, with `method="krylov"`, adaptive Arnoldi (far fewer matvecs at large `t‖A‖`: 2.3–4.4× faster with the default `m=30`, up to 12× with `m=15`, at `t‖A‖ = 4·10³ … 4·10⁵` on a stencil, ADR-0211). Drifted Fokker–Planck and inventory-ladder generators, which `SymmetricOperator` rejects |
| `SpdResolvent` | built by `SymmetricOperator.resolvent(...)` | `.solve(b) -> NDArray`, `.solve_batched(b) -> NDArray` (`b` of shape `(n, nc)`), `.solve_info(b) -> (x, iterations, rel_residual)`, `.method`, `.n` | Factor-once / solve-many `x = (λM + A)⁻¹ b` of an SPD operator, `λ ≥ 0`, `M = diag(mass) > 0`. Tridiagonal: exact `O(n)` LDLᵀ; otherwise PCG with IC(0) or Jacobi. The GIL is released during the solve. PCG is `O(√κ·nnz)` per solve; no multigrid (§62, ADR-0202) |
| `MassKOperator` | `.from_k_and_mass(k_op, m_dense)` | `.evolve(t, v, path="chebyshev", tol=1e-10, m_max=18, n_steps=100) -> NDArray`, `.n()` | Consistent-mass operator `Â = R⁻ᵀ K R⁻¹` where `M = RᵀR`; applies `e^{-t M⁻¹ K}` via Krylov (§55.4) |
| `Etdrk4` | `.from_symmetric_op(op, nonlinearity, h, phi="taylor", tol=1e-10)` (`phi="krylov"`: adaptive Arnoldi φ, ADR-0211) | `.step(u) -> NDArray`, `.integrate(u0, n_steps) -> NDArray` | Cox-Matthews ETDRK4 for `u' = -Au + N(u)`; `"allen_cahn"` nonlinearity `N(u) = u − u³`; arbitrary Python callbacks NOT supported (ADR-0189, §58.3) |

**Example — stiff operator: the rational and implicit paths.** The polynomial
Krylov paths need `O(√(λ_max·t))` mat-vecs; `path="rational"` (shift-and-invert
Lanczos, ADR-0206) reaches `tol` in a fixed number of `O(N)` tridiagonal solves
(27 at `tol = 1e-10`) whatever the stiffness. `path="implicit"` (Jacobi-preconditioned
backward-Euler PCG, accuracy `O(t/n_steps)`, about `√κ` CG iterations per
sub-step with `κ ≈ 1 + (t/n_steps)·λ_max`) handles stiffness such as
`λ_max ≈ 4·10⁷` for smooth data. Rough data (e.g. white noise) at this
stiffness needs more CG iterations than the fixed cap of 400 and raises
`SemiflowError` with `kind == "ConvergenceFailed"`.

```python
import numpy as np
import semiflow

# Stiff 1-D Neumann Laplacian (tridiag(-1, 2, -1), boundary rows 1, -1) * 1e7, N = 400.
n, scale = 400, 1e7
indptr, indices, data = [0], [], []
for i in range(n):
    for j in (i - 1, i, i + 1):
        if 0 <= j < n:
            indices.append(j)
            data.append(scale * ((1.0 if i in (0, n - 1) else 2.0) if j == i else -1.0))
    indptr.append(len(indices))

op = semiflow.SymmetricOperator.from_csr(
    np.array(indptr, dtype=np.int64), np.array(indices, dtype=np.int32),
    np.array(data), n, 1e-6)
v0 = np.linspace(1.0 / n, 1.0, n).reshape(n, 1)
u = op.evolve_batched(1.0, v0, path="implicit", n_steps=100)

# Everything but the constant (null-space) mode has decayed: u -> mean(v0).
assert np.max(np.abs(u - v0.mean())) < 1e-6

# The rational path: 27 tridiagonal solves, error below 1e-10.
u = op.evolve_batched(1.0, v0, path="rational")
assert np.max(np.abs(u - v0.mean())) < 1e-10
```

`mass_lumped_evolve` accepts the same `path=` keyword (`"rational"` uses `M` directly); `MassKOperator.evolve` the polynomial and implicit ones.

### Free functions

`version()`, `heisenberg_heat_kernel()`, `sample_gridfn2d()`, `edge_weight_grad()`, `edge_weight_grad_batched()`, `graph_expmv_frechet()`, `symmetric_op_expmv_frechet()`, `assemble_conservative_csr_1d()`, `mass_lumped_evolve()`, `phi_action()`, `phi_action_batched()`, `phi_combination()`, `shift1d_coeff_grad()`, `richardson_weights()` — see `__init__.pyi` for their signatures.

---

## Performance

GIL release follows the three-phase `py.detach` pattern (ADR-0031):
acquire → snapshot inputs → detach → Rust compute → reacquire. `Send + Sync`
is verified at compile time with `static_assertions`.

Indicative timings on i7-12700K (1000 nodes, 100 steps, `Heat1D`):

| Metric | Value |
|--------|-------|
| Throughput (criterion) | ~56.6 ms per call |
| p99.9 latency (HFT loop, N=1536) | 45 ns/tick |
| Memory footprint | 2.8 MB RSS |

`.with_a_array` uses a pure-Rust `Arc<Vec<f64>>` Catmull-Rom interpolant and
never re-acquires the GIL during `evolve`. `.with_a_function` calls the Python
callables about 10 times per node (`a` twice, `a_prime` 7×, `a_double_prime`
once) per τ-keyed step plan (ADR-0215). There is one plan per `evolve` call of
`n ≥ 2` steps, per chunk of `evolve_chunked` (of ≥ 2 steps) and per Richardson
level, not one per step. Single-step calls, and grids whose plan exceeds 2²²
doubles, pay the callbacks on every step; prefer `.with_a_array` for those.

## Type stubs

`__init__.pyi` and the `py.typed` marker ship with every wheel. Static
type checkers (mypy, pyright, pylance) pick them up automatically.
The `pyrightconfig.json` at the repo root adds `crates/semiflow-py/python`
to `extraPaths` so local development also resolves the stubs correctly
(0 `reportAttributeAccessIssue` errors).
## Bindings

| Language | Package | Notes |
|----------|---------|-------|
| Rust | [`semiflow`](https://crates.io/crates/semiflow) ([docs.rs](https://docs.rs/semiflow)) | The full engine catalogue; `no_std + alloc` |
| Python | [`semiflow-pde`](https://pypi.org/project/semiflow-pde/) | 90 classes and the functions `version()`, `heisenberg_heat_kernel()`, `sample_gridfn2d()`, `edge_weight_grad()`, `edge_weight_grad_batched()`, `graph_expmv_frechet()`, `symmetric_op_expmv_frechet()`, `assemble_conservative_csr_1d()`, `mass_lumped_evolve()`, `phi_action()`, `phi_action_batched()`, `phi_combination()`, `shift1d_coeff_grad()`, `richardson_weights()`; NumPy in/out; complete `.pyi` stubs; see the [PyPI page](https://pypi.org/project/semiflow-pde/) |
| JavaScript / WASM | [`@semiflow/wasm`](https://www.npmjs.com/package/@semiflow/wasm) | Lite build on npm; the heavy-grid engines need a `--features full` build; see the [npm page](https://www.npmjs.com/package/@semiflow/wasm) |
| C / C++ | `semiflow-ffi` | `extern "C"` ABI with `catch_unwind` on every entry point; header `semiflow.h` |

The bindings mirror the user-facing engines, not the internal composition types
(`AxisLift`, `StrangSplit`, … are Rust-side building blocks). Variable
coefficients cross the language boundary as pre-sampled arrays (evaluated in
Rust, no call back into the host language during `evolve`); some classes also
accept a host-language callback (e.g. `Heat1D.with_a_function` in Python,
`Heat1D.withAFunction` in JS, `smf_state_new_with_closure` in C).
The [Bindings guide](https://github.com/VolkovIlia/semiflow/blob/master/docs/BINDINGS.md) has side-by-side examples.

## Documentation

| If you want to… | Read |
|-----------------|------|
| Get started | [Quickstart](https://github.com/VolkovIlia/semiflow/blob/master/docs/QUICKSTART.md) · [User Guide](https://github.com/VolkovIlia/semiflow/blob/master/docs/USER_GUIDE.md) · [Install](https://github.com/VolkovIlia/semiflow/blob/master/docs/INSTALL.md) |
| Use a binding | [Bindings guide](https://github.com/VolkovIlia/semiflow/blob/master/docs/BINDINGS.md) (C / Python / WASM) |
| Browse the Rust API | [docs.rs/semiflow](https://docs.rs/semiflow) |
| See worked examples | [`crates/semiflow/examples`](https://github.com/VolkovIlia/semiflow/blob/master/crates/semiflow/examples/README.md) |
| Understand accuracy and stability policy | [Precision policy](https://github.com/VolkovIlia/semiflow/blob/master/docs/precision-policy.md) · [API stability](https://github.com/VolkovIlia/semiflow/blob/master/docs/api-stability.md) |
| Read design decisions | [Architecture Decision Records](https://github.com/VolkovIlia/semiflow/blob/master/docs/adr/) |
| Follow changes | [CHANGELOG](https://github.com/VolkovIlia/semiflow/blob/master/CHANGELOG.md) · [GitHub Releases](https://github.com/VolkovIlia/semiflow/releases) |

## Accuracy and design

Every numerical claim is gated in CI against a closed-form or high-order
reference oracle (convergence-order and sup-norm tests). Design principles:
`no_std + alloc` core; 3 runtime dependencies; SIMD hot
paths isolated to `src/simd/` (AVX2 on x86_64, NEON on aarch64, scalar
fallback elsewhere); bit-reproducible parallelism across thread counts; no
`unsafe` in the math kernels. Builds with and without `simd` agree to rounding
but not always bit for bit, so golden files are tied to one configuration.

## How to cite

If you use SemiFlow in academic work, cite both the software and the theorem.

```bibtex
@software{volkov2026semiflow,
  author  = {Volkov, Ilia},
  title   = {{SemiFlow}: {Chernoff} Approximation of Operator Semigroups},
  year    = {2026},
  version = {0.15.0-beta},
  doi     = {10.5281/zenodo.20837851},
  url     = {https://doi.org/10.5281/zenodo.20837851}
}

@article{Remizov2025,
  author  = {I. D. Remizov},
  title   = {Chernoff Approximations of the Solution of Linear ODE with Variable Coefficients},
  journal = {Vladikavkaz Math. J.},
  volume  = {27},
  number  = {4},
  pages   = {124--135},
  year    = {2025},
  doi     = {10.46698/a3908-1212-5385-q}
}
```

The concept DOI [10.5281/zenodo.20837851](https://doi.org/10.5281/zenodo.20837851)
always resolves to the latest version; [`CITATION.cff`](https://github.com/VolkovIlia/semiflow/blob/master/CITATION.cff) has the
full entry (GitHub's "Cite this repository" reads it).

## License

Dual-licensed under [MIT](https://github.com/VolkovIlia/semiflow/blob/master/LICENSE-MIT) or [Apache-2.0](https://github.com/VolkovIlia/semiflow/blob/master/LICENSE-APACHE), at your
option.

## Author, contributing, security

SemiFlow is created and maintained by **Ilia Volkov**. See
[CONTRIBUTING.md](https://github.com/VolkovIlia/semiflow/blob/master/CONTRIBUTING.md) and [CODE_OF_CONDUCT.md](https://github.com/VolkovIlia/semiflow/blob/master/CODE_OF_CONDUCT.md).
Report vulnerabilities privately as described in [SECURITY.md](https://github.com/VolkovIlia/semiflow/blob/master/SECURITY.md) —
not in public issues.

