Metadata-Version: 2.4
Name: fast-numerical-integration
Version: 0.2.0
Summary: Fast Cython RK23/RK45/DOP853 integrators that reproduce scipy.integrate.solve_ivp step for step, with a GIL-free C API
Author: Ian Mackey
License-Expression: MIT
Project-URL: Homepage, https://github.com/idm23/fast-numerical-integration
Project-URL: Repository, https://github.com/idm23/fast-numerical-integration
Project-URL: Issues, https://github.com/idm23/fast-numerical-integration/issues
Keywords: ode,numerical-integration,runge-kutta,dop853,cython,solve_ivp,scipy
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: MacOS
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Cython
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
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: Topic :: Scientific/Engineering :: Mathematics
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=2.0
Requires-Dist: scipy>=1.13
Provides-Extra: plot
Requires-Dist: matplotlib>=3.8; extra == "plot"
Dynamic: license-file

# fast-numerical-integration

A bridge for **fast numerical integration** of differential equations in Python.

The core idea: author a system of ODEs in a computer algebra system (SymPy or
similar), **code-print that system to C**, compile it and **bridge it into Python
with Cython**, then integrate it with hand-written numerical integrators
(DOP853, Runge-Kutta families, ...). The goal is to be **considerably faster than
`scipy.integrate.solve_ivp`** while producing reusable, composable tooling that
other projects can drop in.

## Why

`solve_ivp` is flexible and battle-tested, but every right-hand-side evaluation
crosses back into Python, which dominates the cost for stiff or long-horizon
problems. By generating a compiled C RHS and driving it from a Cython integration
loop, we keep the hot path entirely in native code — no per-step Python overhead —
while still authoring the physics symbolically.

## Architecture (target)

```
  SymPy / CAS            Code printing            Cython bridge            Integrators
 ┌────────────┐   ->    ┌──────────────┐   ->   ┌───────────────┐   ->   ┌──────────────┐
 │  dy/dt =   │         │  generated   │        │  compiled .so │        │ DOP853 / RK  │
 │  f(t, y)   │         │  C kernel    │        │  + Python API │        │  step loop   │
 └────────────┘         └──────────────┘        └───────────────┘        └──────────────┘
```

Each stage is meant to be reusable on its own: the code printer, the Cython
build glue, and the integrators are independent tools.

## Current status

**Milestone 1 — SciPy baseline (done).** A `main` entry point solves the classic
[figure-eight three-body choreography](https://en.wikipedia.org/wiki/Three-body_problem#Special-case_solutions)
with `scipy.integrate.solve_ivp` (DOP853) and saves a static trajectory plot.
This gives us a **correct reference solution** and a **timing target** for the
custom engine that follows.

**Milestone 3 — Cython integrators (done).** `RK23`, `RK45`, and `DOP853` are
reimplemented in Cython as a step-for-step copy of SciPy's `RungeKutta`
solvers. They use the same tableaus, initial-step heuristic, error norms, and
step-size controller, so they accept the same steps and make the same number of
RHS evaluations as `solve_ivp`. State-sized work is done in typed-memoryview loops
over buffers that cnumpy allocates once per solve. The step loop makes no NumPy
calls and creates no temporary arrays. Right-hand sides are `cdef` subclasses of
`OdeSystem`, and the integrator calls them through the vtable. A Python callable
also works through an adapter, which is slower.

**Milestone 3b — the rest of `solve_ivp` (done).** `t_eval`, `dense_output`
(an `OdeSolution`), `events` (`terminal` as a bool or a count, plus
`direction`), `vectorized`, `args`, per-component `rtol`/`atol`, `first_step`,
`max_step`, and SciPy's validation and warnings are all supported. Output
points, interpolants, and event times and states match SciPy to round-off,
with identical RHS evaluation counts. That includes DOP853's 3 extra
dense-output stages, which are evaluated only when an interpolant is needed.
Event roots are found with a C port of SciPy's `brentq`. Events can be native
`EventFunction` subclasses, such as the built-in `ComponentCrossing`, so a fully
compiled problem never calls into Python. The stepper classes mirror
`scipy.integrate.OdeSolver` (`step()`, `dense_output()`, `t_old`, `step_size`).
Not supported: complex-valued `y0`.

Package layout (`src/` layout):

```
src/fast_numerical_integration/
  integrators.pyx/.pxd  # GIL-free C core (fni_*), OdeSystem/EventFunction, steppers, solve_ivp, integrate
  batch.pyx             # integrate_batch: many segments in one nogil loop (C API example)
  nbody.pyx             # NBodyGravity: native N-body RHS (stand-in for generated C)
  tableaus.py           # Butcher tableaus, bit-identical to SciPy's
  three_body.py         # figure-eight ICs, three_body_rhs, ThreeBodySystem(.native_rhs())
  visualize.py          # plot_trajectories(solution, out_path) -> PNG
  main.py               # main(): solve with solve_ivp, print stats, write outputs/figure_eight.png
  benchmark.py          # timing harness: scipy vs cython (Python RHS) vs cython (native RHS)
  suite.py              # fni-suite: benchmark + figures for every catalog problem and event scenario
  testing.py            # run SciPy with our summation order (bit-identical comparisons)
  problems/
    catalog.py          # 41 canonical test problems: RHS, exact/published references, invariants
    native.pyx          # native (Cython) RHS and event classes for every catalog problem
    events.py           # event scenarios with exact event times (apsides, pendulum, bouncing ball)
setup.py                # Cython extension build (metadata stays in pyproject.toml)
```

**Milestone 3c — canonical test suite (done).** A catalog of 41 standard
problems, each with a NumPy RHS and an equivalent native one:

- **DETEST** A1–E5, with the published reference end values (Hull et al. 1972;
  Enright & Pryce 1987).
- **Hairer–Nørsett–Wanner:** Arenstorf, Pleiades, Brusselator.
- **Long-time:** Kepler over 100–1000 orbits, the outer solar system over
  2×10⁵ days, the figure-eight choreography over 50 periods.
- **Stability:** Prothero–Robinson, Van der Pol for μ = 1…1000, Robertson,
  Lorenz, finite-time blow-up.
- **Event scenarios:** Kepler apsides, a large-amplitude pendulum, a bouncing
  ball with restarts.

The tests are grouped by the oracle they check against:

| directory | checks |
|---|---|
| `tests/parity/` | SciPy's own `solve_ivp` tests, ported; RHS-evaluation parity on every catalog problem; **bit-identical** trajectories against SciPy patched to our summation order |
| `tests/accuracy/` | DETEST global error vs tolerance; observed order 3/5/8 with fixed steps; dense-output order; references cross-checked against Radau and LSODA |
| `tests/long_time/` | Kepler: linear invariant drift and quadratic phase-error growth; two independent solar-system formulations agree; linear invariants exact to round-off |
| `tests/stability/` | step size rides the stability boundary computed from each tableau (2.51 / 3.31 / 6.39); stiffness scaling; blow-up; Lorenz error growth at the Lyapunov rate |
| `tests/events/` | exact event times on physical problems; edge cases (t0, t_end, simultaneous/terminal ordering, root pairs and grazing touches missed exactly as in SciPy) |
| `tests/baseline/` | SciPy baseline, benchmark harness, suite smoke test |

Findings worth knowing:

- RK23 and RK45 match SciPy's RHS-evaluation count on every problem. DOP853
  can differ by a few percent: its error estimate is a heavily cancelling
  sum, so near round-off level summation order flips accept/reject decisions.
  SciPy shows the same sensitivity to a 1-ulp change in `y0`. With SciPy
  patched to our summation order (`fast_numerical_integration.testing`),
  every trajectory is bit-identical.
- RK45 at rtol=1e-3 goes unstable on Robertson (|y| ~ 1e7, then fails). This
  is a known Dormand–Prince stiffness pathology, and SciPy does the same.

### Using the Cython engine

`integrators.solve_ivp` takes the same arguments as `scipy.integrate.solve_ivp`
(default method `RK45`):

```python
import numpy as np
from fast_numerical_integration.integrators import ComponentCrossing, solve_ivp
from fast_numerical_integration.three_body import figure_eight_system

system = figure_eight_system()
sol = solve_ivp(
    system.native_rhs(), (0.0, system.period), system.y0, method="DOP853",
    t_eval=np.linspace(0.0, system.period, 2000), dense_output=True,
    events=[ComponentCrossing(1, direction=1)],    # y1 crosses 0 upward
    rtol=1e-10, atol=1e-12,
)
sol.t, sol.y              # solution at t_eval
sol.sol(1.234)            # dense interpolant anywhere in the interval
sol.t_events, sol.y_events
```

Plain Python RHS and event callables also work, with SciPy's `terminal` and
`direction` attributes. They are slower because every call goes back through
Python.

### Native systems, runtime parameters and the C API

A native RHS is a GIL-free C function plus a context pointer. To write one,
subclass `OdeSystem` (or `EventFunction` for events) in a `.pyx` file. In
`__init__`, register the C function with `c_fn` / `c_ctx` and declare the
runtime parameters:

```cython
from fast_numerical_integration.integrators cimport OdeSystem

cdef int _vdp(double t, const double* y, double* dydt, Py_ssize_t n, void* ctx) noexcept nogil:
    cdef const double* p = <const double*> ctx          # p[0] = mu
    dydt[0] = y[1]
    dydt[1] = p[0] * (1.0 - y[0] * y[0]) * y[1] - y[0]
    return 0

cdef class VanDerPol(OdeSystem):
    def __init__(self, double mu=1.0):
        self.n_states = 2
        self._init_params([mu], ("mu",))
        self.c_fn = _vdp
        self.c_ctx = self._params_ptr()
```

Parameters change between solves without recompiling, which is what
continuation needs:

```python
vdp = VanDerPol()
for mu in (1.0, 10.0, 100.0):
    vdp.params = [mu]
    sol = solve_ivp(vdp, (0, 20), [2.0, 0.0], method="DOP853")
```

(Overriding `cdef eval` still works, but such systems need the GIL.)

All three entry points run on one C core, so they take bit-identical steps:

| entry point | what it does |
|---|---|
| `solve_ivp(...)` | full SciPy-compatible API: trajectory, `t_eval`, dense output, events |
| `integrate(...)` | endpoint + event locations only (shooting). Releases the GIL when the system and events are native, so Python threads run in parallel |
| `fni_integrate(...)` (C) | `cimport`able from other Cython modules and callable inside `with nogil`. `batch.integrate_batch` is the reference example: it propagates many segments at ~1 µs per segment |

```cython
from fast_numerical_integration.integrators cimport (
    fni_integrate, fni_options, fni_stats, FNI_DOP853, FNI_OK)

with nogil:
    rc = fni_integrate(FNI_DOP853, rhs_fn, ctx, n, t0, y0, t1, &opt,
                       events, n_events, y_out, rec_event, rec_t, rec_y, max_records, &stats)
```

The lower-level `fni_solver_init` / `fni_solver_step` / `fni_solver_dense_eval` /
`fni_events_check` functions expose single steps for custom drivers. All
solver state lives in a caller-owned struct, so concurrent integrations are
safe if the RHS kernels are reentrant. Every catalog kernel in
`problems/native.pyx` is written that way.

## Installation

```bash
pip install fast-numerical-integration            # or: uv add fast-numerical-integration
pip install "fast-numerical-integration[plot]"    # + matplotlib, for fni / fni-suite / visualize
```

Pre-built wheels cover Linux x86_64 and macOS on Apple Silicon for CPython
3.11–3.14, and every wheel passes the full test suite before it's
published. Other platforms install from the source distribution, which needs
a C compiler. All extensions are built with floating-point contraction off,
so results are identical across platforms.

To `cimport` the C API from your own Cython package, add
`fast-numerical-integration` and `cython` to your `[build-system] requires`.
`examples/downstream/` is a minimal working package; the release workflow
builds it against every release.

## Quickstart (development)

Requires [`uv`](https://docs.astral.sh/uv/).

```bash
uv sync                                        # create the venv, install deps, compile Cython
uv run fni                                     # or: uv run python -m fast_numerical_integration.main
```

This writes `outputs/figure_eight.png` and prints solver stats, e.g.:

```
System:        figure-eight
Method:        DOP853 (scipy solve_ivp)
Steps taken:   2000 sampled, 1454 RHS evals
Wall time:     22.30 ms
Energy:        E0=-1.287142e+00  ET=-1.287142e+00  |drift|=2.61e-10
Saved plot ->  outputs/figure_eight.png
```

The tiny energy drift over one period confirms the dynamics and initial
conditions are correct.

The Cython extensions build when the package installs. `uv` tracks the `.pyx` and
`.pxd` files as cache keys, so `uv sync` or `uv run` rebuilds them after an
edit. Building needs a C compiler.

### Benchmark

```bash
uv run fni-bench                               # or: uv run python -m fast_numerical_integration.benchmark
```

Times `RK23`, `RK45`, and `DOP853` on the figure-eight system (one period,
rtol=1e-10, atol=1e-12) with three engines:

- `scipy`: `scipy.integrate.solve_ivp`.
- `cython+pyrhs`: the Cython integrators calling the NumPy RHS through Python.
  The RHS cost matches SciPy's, so this measures integrator overhead alone.
- `cython`: the Cython integrators with the native `NBodyGravity` RHS.

Representative run:

```
solver                 best (ms)      mean±sd (ms)    nfev   eval/ms     E drift   period err
---------------------------------------------------------------------------------------------
scipy RK23               626.758    636.353±10.984   43064      68.7    7.66e-10     6.93e-08
cython+pyrhs RK23        432.419     435.549±2.256   43064      99.6    7.67e-10     6.93e-08
cython RK23                3.747       3.966±0.299   43064   11494.1    7.67e-10     6.93e-08
scipy RK45                39.150      39.270±0.093    2936      75.0    2.98e-10     7.78e-08
cython+pyrhs RK45         29.834      31.064±1.231    2936      98.4    2.98e-10     7.78e-08
cython RK45                0.201       0.207±0.005    2936   14633.5    2.98e-10     7.78e-08
scipy DOP853              15.282      15.353±0.060    1202      78.7    2.61e-10     7.29e-08
cython+pyrhs DOP853       12.178      12.234±0.037    1202      98.7    2.61e-10     7.29e-08
cython DOP853              0.086       0.091±0.003    1202   13906.7    2.61e-10     7.29e-08
```

The Cython engine makes exactly the same number of RHS evaluations as SciPy and
reaches the same accuracy. With a native RHS, DOP853 runs about **175x faster**
than `solve_ivp`. With the Python RHS the gain is only ~1.3x, so most of SciPy's
time is spent calling the Python RHS. That is why the remaining milestones
generate compiled kernels. Numbers vary by machine.

### Benchmark and visualize every problem

```bash
uv run fni-suite                    # all problems + event scenarios, ~2 minutes
uv run fni-suite --problems D5 pleiades --quick
```

This writes one figure per problem to `outputs/suite/<group>/<name>.png`. Each
figure shows the solution, error along the trajectory (or conserved-quantity
drift, or step size), a DOP853 work-precision curve and wall-time bars for
SciPy and Cython. It also writes one figure per event scenario and
`summary_speedup.png`, plus `summary.md` / `summary.csv` tables. Cython's
end-to-end speedup with DOP853 ranges from about 6× (C4, a dense 51×51 linear
system, where the work is the same matrix-vector product as NumPy's) to about
500× (outer solar system), with identical step counts and accuracy.

Run the tests:

```bash
uv run pytest               # ~1.5 minutes
uv run pytest --runslow     # include long integrations
uv run pytest tests/events  # one oracle group
```

CI (`.github/workflows/tests.yml`) runs the full suite, including
`--runslow`, on Python 3.11, 3.12 and 3.13 for every pull request and push
to `main`. The `all-tests` check aggregates the matrix and is the status to
require before merging.

## Roadmap

- [x] **Milestone 1** — SciPy `solve_ivp` figure-eight baseline + visualization.
- [x] **Milestone 2** — Benchmark harness timing `solve_ivp` across `RK23`,
      `RK45`, and `DOP853` on the same system (baseline numbers to beat).
- [x] **Milestone 3** — Cython RK23 / RK45 / DOP853 with adaptive steps, matching
      SciPy step for step; native `OdeSystem` RHS interface.
- [x] **Milestone 3b** — Dense output, `t_eval`, events, and the rest of
      SciPy's `solve_ivp` options for these methods.
- [x] **Milestone 3c** — Canonical test problems (DETEST, Hairer, long-time,
      stability, events), split by oracle, with a benchmark/figure suite.
- [x] **Milestone 3d** — GIL-free C core with a `cimport`able C API
      (`fni_integrate`, single-step functions), runtime `params` on native
      systems and events, `integrate()` endpoint API, `integrate_batch`.
- [ ] **Milestone 4** — SymPy → C code printing for the RHS.
- [ ] **Milestone 5** — Build glue that compiles the generated C into an
      `OdeSystem` subclass; benchmark against the baseline.

## Releasing

Releases are published from GitHub Actions (`.github/workflows/release.yml`)
with PyPI trusted publishing, so no tokens are stored.

1. Bump `version` in `pyproject.toml` and merge the change to `main`.
2. Optional rehearsal: run the *Release* workflow manually with
   `target: testpypi`, then `pip install -i https://test.pypi.org/simple/
   --extra-index-url https://pypi.org/simple/ fast-numerical-integration`.
3. Tag and push: `git tag vX.Y.Z && git push origin vX.Y.Z`. The workflow
   checks that the tag matches the version, builds and tests every wheel, and
   uploads to PyPI.

## License

MIT. See [LICENSE](LICENSE).
