This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!
viv-lang.org
Python Viv

vp.ode.batch

Solve one initial value problem for many parameter sets and initial states in one call.

Part of the distribution vipyr; it runs from a checkout of the repository. Planned on PyPI as pip install vipyr.

The case SciPy's solve_ivp handles worst: thousands of small systems, one per parameter set (a dose-response sweep, a posterior sample, an ensemble of epidemics). A Python loop pays solve_ivp's per-call and per-step overhead for each; vp.ode.batch runs the compiled solver over all the rows in one call. Each row is an independent adaptive solve with its own steps, reported on the common grid t_eval.

Signature

vp.ode.batch(rhs, t_span, y0, params, t_eval, *, method='auto', rtol=0.001, atol=1e-06, events=None, direction=None, terminal=None, jac=None, first_step=None, max_step=None, max_steps=1000000, fast=False)

Parameters

rhs : str | function, required
The right-hand side f(t, y, p), as for vp.ode.solve.
t_span : (float, float), required
The interval (t0, t1), the same for every solve.
y0 : ndarray[float64], required
The initial state, length n, or one row per solve (B, n).
params : ndarray[float64], required
One parameter vector per solve, shape (B, k): row b is the p of solve b.
t_eval : ndarray[float64], required
The output times, the same for every solve.
method : str, default 'auto'
As for vp.ode.solve.
rtol : float, default 0.001
Relative tolerance.
atol : float | ndarray[float64], default 1e-06
Absolute tolerance, one number or one per component.
events : str | function, default None
As for vp.ode.solve; each solve reports its first event.
direction : ndarray[float64], default None
As for vp.ode.solve.
terminal : ndarray[bool], default None
As for vp.ode.solve; after a terminal event the outputs are NaN.
jac : str, default None
As for vp.ode.solve.
first_step : float, default None
As for vp.ode.solve.
max_step : float, default None
As for vp.ode.solve.
max_steps : int, default 1000000
As for vp.ode.solve, per solve.
fast : bool, default False
As for vp.ode.solve.

Returns

A frozen OdeBatch. Iterating it gives the fields below in order; it also has params, tests, log and check().

t : ndarray[float64]
The output times (t_eval).
y : ndarray[float64]
The solutions, shape (B, n, len(t)).
status : ndarray[int64]
Each solve's status, as for vp.ode.solve.
success : bool
Every solve reached t1 or a terminal event.
first_event : ndarray[float64]
Each solve's first event time, NaN if none (None without events).
stats : dict
Per solve: nfev, njev, nlu, steps, rejected (arrays of length B).

Raises

ValueError
A method, interval, tolerance or array shape the solver cannot take; vipyr._ode.OdeProblemError (a ValueError) when the rhs, events or jac cannot be made into a Viv problem, or py2viv refuses a Python rhs.
vipyr.VivCompileError
A SyntaxError: the rhs's .viv file does not compile.
TypeError
An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError
An ImportError: no C compiler, or the program's C does not build.

See also

  • vp.ode.solve : Solve an initial value problem y' = f(t, y, p), stiff or not, with the right-hand side in Viv.

Notes

Row b of the result is exactly what vp.ode.solve(rhs, t_span, y0[b], params=params[b], t_eval=t_eval, ...) returns, bit for bit: the batch is a loop in Viv over the same solver, not a vectorized approximation.

The output is y[b, i, j]: solve b, component i, time t_eval[j]. A solve stopped by a terminal event has NaN after it; first_event[b] is its first event time. status[b] and the per-row stats say how each solve ended and what it cost.

References

  1. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer (1996). doi:10.1007/978-3-642-05221-7
  2. H. Bateman, The solution of a system of differential equations occurring in the theory of radio-active transformations, Proc. Cambridge Philos. Soc. 15, 423-427 (1910): the closed form the pharmacokinetic batch is checked against.

Examples

import numpy as np
import vipyr as vp
# 1000 SIR epidemics, R0 from 1.2 to 4, recovery in 7 days
beta = np.linspace(1.2, 4.0, 1000) / 7.0
P = np.column_stack([beta, np.full(1000, 1.0 / 7.0)])
r = vp.ode.batch("sir", (0.0, 180.0), [0.999, 0.001, 0.0], P, t_eval=np.linspace(0, 180, 181),
                 method="dopri5", rtol=1e-8, atol=1e-10)
print(r.y.shape, round(r.y[-1, 2, -1], 4))     # everyone who was ever infected, at R0 = 4
assert r.check()
(1000, 3, 181) 0.9802

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference
40 one-compartment pharmacokinetic models against Bateman's closed form, worst error 2.8 × 10−10 mg/L at rtol 10−9 (lib/ode/nonstiff_events.viv); every row equal to vp.ode.solve's to the bit (python/tests/test_ode.py).

Source