API reference

Generated from the docstrings of the installed package. For a short list of the functions most specs need, see the API overview.

Core: checks, results, order verdicts

physverify.checks

Result types shared by every check.

Result dataclass

failed property

Neither a pass nor a skip: FAIL, ERROR or INCONCLUSIVE.

Check dataclass

One runnable check. family says which kind of evidence it gives:

mms manufactured-solution order study against an independent reference operator reference comparison with an independently written reference solution invariant drift of a conserved quantity as the resolution is refined limit a known limiting case (analytic range, normal-mode frequencies, ...) patch a manufactured solution the scheme should reproduce exactly (to round-off), so any error at all is a bug, whatever its order

SkipCheck

Bases: Exception

Raise inside a check to report SKIP with this message (e.g. an optional dependency is missing).

requires(*modules)

Raise SkipCheck unless every module imports, so a missing optional dependency is a SKIP, not a pass.

within_budget(resolutions)

Yield the resolutions in order, stopping early when the next one would not finish inside the time budget.

The cost of the next run is guessed as the last run's time scaled by the ratio of resolutions. The order study then fits what it has; fewer than two usable resolutions give INCONCLUSIVE. Outside a budget (no budget_s) every resolution is yielded.

run_check(check, budget_s=None)

Run one check. With a budget (check.budget_s, else budget_s), order studies stop adding resolutions when the next would overrun it (see within_budget), and a check still running at twice the budget is reported INCONCLUSIVE so the later checks still run. Python cannot stop a thread, so such a check keeps running in the background until the process exits; its result is discarded.

physverify.order

Observed order of accuracy from errors at several resolutions.

pairwise_orders(hs, errors)

p_i = log(e_i / e_{i+1}) / log(h_i / h_{i+1}) for consecutive resolutions.

fit_order(hs, errors)

Least-squares slope of log(error) against log(h).

fit_order_error(hs, errors)

(order, standard error): the least-squares slope of log(error) against log(h) and its standard error from the scatter about the line (needs at least three points; two give an error of nan). A fitted 3.9 with standard error 0.05 supports order 4; 3.9 with 0.6 supports nothing yet.

gci(f_fine, f_medium, f_coarse, r, *, safety=1.25)

Grid convergence index of three solutions on grids refined by the constant ratio r (Roache 1994; Celik et al., J. Fluids Eng. 130, 078001, 2008): the observed order p = ln(e32/e21) / ln r, the Richardson extrapolate, and GCI_fine = safety |e21 / f_fine| / (r^p - 1), a relative error band on the fine solution. Oscillatory convergence (e21 and e32 of opposite sign) is reported, and then p is nan. So is the case where the change between grids does not shrink on refinement (not_converging): the formula would give p <= 0 and a negative or huge error band, which means nothing.

effective_floor(floor, scale, rel=1e-13)

Round-off grows with the size of the solution, so the noise floor does too.

The floor is the larger of the absolute floor (right for solutions of size up to about 10) and rel times the solution's size (about 500 machine epsilons, room for round-off accumulated over a few hundred steps). A non-finite scale (a diverged run) leaves the floor unchanged so that the non-finite error is still reported.

check_exact(hs, errors, *, floor=1e-12, scale=1.0)

Require the error to be at round-off level at every resolution.

For a manufactured solution that lies inside the scheme's discrete solution space (a "patch test": for example, polynomials a finite-difference stencil differentiates exactly), a correct code reproduces it to round-off whatever the step sizes. Order of convergence cannot see an error that is as small as the scheme's own truncation error and shrinks at the same rate; this can, because the expected error is zero rather than small.

EXACT every error is at most floor in size FAIL some error is larger, or non-finite

scale is the size of the solution; see effective_floor.

check_order(hs, errors, expected, *, tol_low=0.3, tol_high=1.0, floor=1e-12, fit_points=4, scale=1.0, tol=None, atol=None)

Compare the observed order of convergence with the order the code claims.

PASS fitted order lies in [expected - tol_low, expected + tol_high] EXACT every error is below floor, so the order cannot be measured FAIL order too low or too high, error plateau, or non-finite error INCONCLUSIVE fewer than two resolutions are above the noise floor

Errors may be signed (computed minus exact): their sizes are used, so a large negative error is not taken for round-off.

floor is the error level below which round-off dominates, for a solution of size about 1. Pass scale (the size of the solution, for example the largest absolute value of the exact answer) and the floor grows with it.

tol sets tol_low and tol_high to the same value (a common way to write it; without the alias an adapter that wrote tol= failed with a TypeError after its whole run).

atol adds an accuracy requirement on the finest resolution. It is reported separately: when the order is as claimed and only atol is missed, the message says so, because that is a question of how fine the grid must be, not a convergence bug (batch 7 had an adapter whose absolute tolerance failed while the observed order was exactly as claimed).

Order studies and manufactured solutions

physverify.studies

Order studies other than MMS, and a helper that bundles them for an ODE code.

reference_study(integrate, rhs_ref, y0, t_end, ns, expected_order, ref_factor=16, **kw)

Compare with a reference solution from an independent right-hand side.

The reference is RK4 on rhs_ref with ref_factor times the finest step count, so its own error is negligible next to the code under test.

invariant_study(integrate, quantity, y0, t_end, ns, expected_order, **kw)

Relative drift of a conserved quantity at t_end as the step shrinks.

Drift is only required to shrink at least as fast as the scheme's order. RK4 on a pendulum, for example, drifts in energy at about order 5, so a drift that converges faster than expected is not a failure here.

linear_manufactured(rhs_ref, y0, slope)

u(t) = y0 + slope * t, with its forcing taken from the reference operator.

patch_study(integrate, rhs_ref, y0, slope, t_end, ns, *, floor=1e-11)

Exactness test for an ODE code: a manufactured solution linear in t.

Every consistent explicit Runge-Kutta method (Euler, midpoint, RK4, ...) reproduces u(t) = y0 + slope * t exactly, whatever the step: each stage lands on u at its own stage time, so each stage derivative is exactly slope, and weights that sum to one return y0 + slope * (t + h). So any error above round-off is a bug, including bugs whose error shrinks at the scheme's own order (an extra step for Euler, a wrong stage time) that every order-of-convergence check accepts. See docs/patch.md.

It does not see a wrong but consistent stepper (Euler in place of RK4 is still exact here): the order studies cover that.

smoke_check(label, integrate, rhs_ref, y0, t_end, n, *, rtol=0.001)

What an ordinary unit test would do: one run, compared with a reference within a tolerance.

Not one of the verification families. It is here so that the benchmark can report what a plain test would have caught (see the benchmark comparison, pvbench.compare).

ode_checks(label, integrate, *, rhs_ref, mms, y0_free, t_end, ns, order, quantity=None, floor=1e-12, patch=None, integrate_trajectory=None)

The standard bundle for an ODE code: MMS, reference solution, invariant.

patch (y0, slope) of a linear manufactured solution; adds the exactness check integrate_trajectory same contract as integrate but returning every state; adds an MMS check on the largest error over the whole run

time_grid_check(times_of_run, t_start, t_stop, ns, *, tol=1e-09)

The run stops at t_stop after exactly n steps, for step sizes that divide the interval.

times_of_run(n) runs the code with step (t_stop - t_start) / n and returns the times it reports. A loop such as while t <= t_end: t += h takes one step too many whenever t lands exactly on t_end in floating point (h = 0.5, 0.25) and the right number otherwise (h = 0.1), so include step sizes that are exact in binary. An order study that compares at the time nearest t_stop cannot see this; this check can.

backward_mms_study(integrate_back, mms, t_end, ns, expected_order, **kw)

MMS order study integrating from t_end back to 0.

integrate_back(source, y_end, t_end, n) integrates y' = f(t, y) + source(t) from t_end down to 0 in n steps and returns y(0). Use a manufactured solution whose forcing depends on t: a code that evaluates its stages at t + h/2 while stepping backwards still passes on autonomous problems but drops to first order here.

balance_check(measure, expected, resolutions, *, tol=1e-09, label='')

Check an exact discrete balance at every resolution: measure(n) must equal expected to tol.

A conservative scheme keeps a budget exactly on every grid: the total heat with insulated walls, the heat gained from a known boundary flux or source over a run. Three real codes failed such a budget while still converging: a zero-flux wall made by copying the neighbour value lost 8.5% of the heat, a source added without the factor dt depended on the number of steps, and boundary nodes given full cells took in 7% too much heat. A convergence study alone can pass these, so the check reports whether the mismatch shrinks with resolution.

size_sweep_check(run_at_size, sizes=(15, 16, 31, 32, 63, 64), *, label='')

Call run_at_size(n) for odd and even sizes; FAIL if it raises for some of them.

Spectral and FFT-based codes often build their wavenumber arrays by splitting the grid in two halves, which only works for even n. A real code raised a ValueError for every odd grid size while being exact for even ones; an order study on 16, 32, 64, 128 points cannot see that. run_at_size may also return False to mark a wrong result at that size.

damped_sinusoid_fit(t, y, t_start, t_stop)

Complex frequency omega_R - i gamma of y ~ A exp(-gamma t) cos(omega_R t + phi) on [t_start, t_stop].

Two-pole Prony fit (linear prediction y_n = a1 y_{n-1} + a2 y_{n-2} by least squares on a uniform grid), so no starting guess is needed. Start the window after the direct pulse has passed the observer and stop before the late-time tail or reflections dominate; a window that includes the pulse gave a wrong frequency for a real ringdown code (J216). Returns omega_R - i gamma with omega_R >= 0, to compare with tabulated quasinormal modes (references.schwarzschild_qnm) or any damped oscillator.

trajectory_order_study(runs, exact, t_compare, expected_order, *, period=None, time_tol=1e-09, **kw)

Order study from whole trajectories, compared at one fixed time, guarding against the usual study-design pitfalls (batch 7 of the corpus study).

runs is a list of (t, y) pairs, one per step size, as the code under test returns them; exact(t) is the solution. Each run is compared at t_compare, found on its own time grid; the step size is the run's mean step.

Pitfalls it catches: * A run whose grid does not contain t_compare: a loop over int(t_end / h) steps stops one step early (int(2 pi / (2 pi / 200)) = 199), which made an order-4 method look first order (J348). FAIL, naming the run and where it ended; it is the study, or the code's step count, that is wrong. * Comparing arrays index by index, or a maximum over the last N samples, moves the comparison time with h (J155). Comparing at one time on every grid avoids both. * t_compare a whole number of periods: errors in the phase can cancel there, and a symplectic Euler showed order 2 (J155). Noted in the message, with the advice to also compare at a time that is not a multiple.

coverage_check(simulate, interval, prior_sampler, *, level=0.68, trials=400, seed=0, alpha=0.001, label='')

Frequentist coverage of an inference code: the stated level interval should contain the truth in that fraction of repeated experiments.

Each trial draws the truth from prior_sampler(rng), simulates data with simulate(theta, rng) and asks interval(data, level) for (low, high). The truth must come from the same prior the analysis assumes: drawing it from a narrower range gave a spurious coverage of 0.865 for a nominal 0.68 interval (J399), which is why the prior sampler is a required argument and not a fixed truth. For a frequentist (not Bayesian) interval any fixed truth works; pass a sampler that returns it.

FAIL when an exact binomial test rejects the stated level at alpha (two-sided; the false-alarm rate of a correct code is at most alpha). Scalar parameters; for vectors, check each component separately.

decay_rate(t, y, *, t_skip=None, windows=8, steady_tol=0.05)

Exponential decay rate gamma of |y| ~ A exp(-gamma t), and the time the fit starts.

A run started from a state that is not a single decaying mode has a start-up transient (a shear wave started from equilibrium, a heat profile with higher harmonics); fitting from t = 0 gives the wrong rate (batch 7). With t_skip the fit uses t >= t_skip. Without it, the record is cut into windows pieces, the local rate is fitted on each, and the fit starts at the first piece after which every local rate stays within steady_tol (relative) of the last one. Returns (gamma, t_start). What is left of the transient after the start biases the rate by up to about steady_tol; lower it (and use more windows) for a tighter fit, or pass t_skip.

y should be a positive amplitude (an energy, a mode amplitude, an envelope); for an oscillating signal use damped_sinusoid_fit. Raises ValueError when no steady stretch is found.

residual_check(solve, residual, tols=(0.0001, 1e-06, 1e-08), *, slack=10.0, label='')

An iterative solver that is asked for tolerance tol returns an answer whose residual is about tol.

solve(tol) returns the solution, residual(x) its residual computed independently (for example the norm of b - A x). PASS when residual <= slack * tol for every tol. A residual that does not follow the requested tolerance is the sign of a solver that does almost no work: a real pressure solver did one Jacobi sweep whatever the tolerance, and its answers looked plausible.

physverify.mms

Method of Manufactured Solutions for systems of ODEs.

You pick a smooth function u(t). Substituting it into an independent reference operator f_ref leaves a residual s(t) = u'(t) - f_ref(t, u(t)). The code under test is then run on y' = f_test(t, y) + s(t). If f_test agrees with f_ref, u(t) is its exact solution and the error must shrink at the scheme's design order. If f_test differs, the error stops shrinking.

Two points that are easy to get wrong:

• The reference operator must be written independently of the code under test. If s(t) is built from the code's own right-hand side, a wrong right-hand side passes.
• The error is measured at a fixed final time, so the code must be driven for an exact number of steps; see the integrate contract below.

The integrate contract used throughout:

integrate(source, y0, t_end, n_steps) -> state at t_end

source(t) returns the additive forcing (a vector, or 0.0). The adapter for a given code decides how to inject it.

complex_step_derivative(u, t, h=1e-30)

du/dt to machine precision, for u written with NumPy functions of t.

u must accept a complex t and avoid abs(), max() and similar non-analytic operations.

mms_order_study(integrate, mms, t_end, ns, expected_order, **order_kwargs)

Run integrate with the manufactured forcing at several step counts.

mms_trajectory_study(integrate_trajectory, mms, t_end, ns, expected_order, **order_kwargs)

Like mms_order_study, but the error is the largest one over the whole run.

integrate_trajectory(source, y0, t_end, n_steps) returns the n_steps + 1 states at t_k = k * t_end / n_steps, starting with y0. A code whose error happens to cancel at the final time, or which goes wrong part way through and recovers, is caught here and not by the final-state check.

physverify.symbolic

Manufactured solutions with symbolic derivatives (optional; needs SymPy).

pip install physverify[symbolic]

Hand-derived source terms are a common place for the verification itself to go wrong. Here you write the manufactured solution and the PDE operator once, symbolically, and the source term is derived and turned into NumPy functions.

Example, the 1D heat equation u_t = alpha u_xx + s:

import sympy as sp
x, t = sp.symbols("x t")
u = sp.exp(-t) * sp.sin(sp.pi * x)
ms = manufactured_pde(u, lambda u: sp.diff(u, t) - 0.7 * sp.diff(u, x, 2), (x, t))
ms.source(xs, 0.3)   # NumPy array of s at the points xs, time 0.3

The operator must be written independently of the code under test, for the same reason as the ODE reference operator (see physverify.mms).

manufactured_pde(u_expr, operator, variables)

Source s = operator(u) for a manufactured solution u.

operator(u) returns the PDE residual for u, written with sympy.diff, so that u solves operator(u) = s exactly. variables is the tuple of SymPy symbols the NumPy functions take, in order (for example (x, t) or (x, y, t)).

manufactured_ode(u_exprs, rhs_ref, t_symbol)

A physverify.ManufacturedODE whose derivative is taken symbolically.

u_exprs is a list of SymPy expressions in t_symbol, one per state component.

Identities, stencils, fluxes, operators

physverify.identities

Checks that need no exact solution: round trips, derivatives, sibling agreement and identities.

These are the checks that found most of the bugs in the batch-8 corpus study (popular libraries): two functions that should agree but do not (an air-absorption term with the opposite sign in one of two reverberation formulas), an operation that should be undone by its inverse but is not (a zero-time rotation that shifts an image by half a pixel), an analytic derivative compared with finite differences, and textbook identities (Poisson's equation, unitarity, conservation).

Each function returns a Result (PASS or FAIL with the worst case in the message); wrap it in a Check("limit", label, lambda: ...) to run it with the others. Values may be real or complex scalars or arrays. The reference side must come from the physics, not from the code under test: a round trip or a sibling comparison can pass when both sides share the same mistake.

check_roundtrip(forward, inverse, inputs, *, rtol=1e-09, atol=0.0)

inverse(forward(x)) == x for every x in inputs.

Catches conversions that disagree with their own inverse: a reference pixel off by half a pixel, a dewpoint formula that does not invert the saturation formula, a unit conversion with the wrong constant on one side. Both functions come from the code under test, so a PASS says only that they are consistent, not that either is right.

check_agree(a, b, inputs, *, rtol=1e-09, atol=0.0, names=('a', 'b'))

Two implementations that must give the same answer do: a(x) == b(x) for each x in inputs.

Use it for sibling code paths (a 2D and a 3D version, a fast path and a slow path, the min and max variants, two formulas that coincide in a limit). Each element of inputs is a tuple of arguments, or a single argument.

fd_derivative(f, x, *, h=None, order=1)

Central finite-difference derivative of f at x (scalar x), with one Richardson step.

order is 1 or 2 (first or second derivative). The default step is cbrt(machine epsilon) times max(1, |x|) for the first derivative and its square root for the second; the Richardson step combines h and h/2 to cancel the leading error term.

check_derivative(f, df, points, *, order=1, h=None, rtol=1e-06, atol=1e-10)

An analytic derivative df agrees with a finite-difference derivative of f at each point.

This is how a thermodynamic-consistency claim is settled cheaply: a reader's derivation said an activity-coefficient derivative dropped a term, and comparing it with finite differences of the activity coefficients themselves showed agreement to 1e-10, so the claim was wrong. The default tolerance allows for the finite-difference error; tighten it only with a well-scaled h.

check_identity(cases, *, rtol=1e-09, atol=0.0, what='identity')

Each (label, lhs, rhs) in cases satisfies lhs == rhs.

For identities written from the physics: Poisson's equation (Laplacian of the potential equals 4 pi G rho), unitarity of a lossless scattering matrix (S^H S = I), E^2 = p^2 + m^2, a limit that must reproduce a simpler formula. Compute the right-hand side yourself, not with the code under test.

laplacian_fd(f, point, *, h=0.001)

Cartesian Laplacian of a scalar field f(x) (x a 1D array) by second-order central differences.

check_unitary(matrix, *, atol=1e-10)

A matrix (or a stack of matrices, last two axes) is unitary: M^H M = I.

A lossless network's scattering matrix, a time-evolution operator, a rotation.

vary_one(test, defaults, changes, *, rtol=1e-09, atol=0.0)

Re-run a test that passes at the default parameters with one parameter changed at a time.

test(**params) returns a Result, a bool, or a (got, ref) pair with the reference from the physics. defaults are the usual values (often symmetric ones: Lx = Ly, nx = ny, zero boundary data, s polarisation, vacuum); changes maps a parameter to a value, or a list of values, to try instead while the others keep their defaults. FAIL names each change that breaks the test.

In the 2026-10 corpus sessions 7 of the confirmed bugs passed every check at the example values and appeared only when one parameter left them: a y-derivative scaled with Lx (only when Lx != Ly), a Neumann row that dropped the boundary value (only when g != 0), the p-polarisation path of a grating solver, an inertia formula right for rectangles only, a factor that cancels for low-loss media.

check_scale(got, ref, *, rtol=1e-06, min_correlation=0.999)

The output has the right size, not only the right shape: got = k * ref with k = 1.

Fits the single factor k (least squares) and the correlation of got with ref. A correlation near 1 with k != 1 is a scale or units error (a missing 1/(2 dx), a missing time step, a factor 1/2); the message names a nearby simple factor (2, pi, sqrt 2, ...) or gives k to match against a grid spacing or a sampling rate. In the 2026-10 corpus sessions the ratio alone pointed to the bug in four repositories. A low correlation means the shape itself is wrong, and a plain comparison says more.

check_vec_convention(superop, f, x, *, expected=None, atol=1e-10)

Which vectorisation a superoperator follows: superop @ vec(X) == vec(f(X)), with vec in C (row-major) or F (column-major) order.

Kronecker constructions of Liouvillians and other superoperators depend on the convention (vec(A X B) = (B^T kron A) vec(X) holds for column stacking); code that builds S for one and reshapes with the other gives a wrong but well-formed answer (corpus repo L204). With expected ('C' or 'F'), FAIL unless S follows it; without, PASS when S follows either, and say which.

physverify.stencils

Finite-difference stencils: consistency and order from the weights alone (no grid, no solver).

A stencil for the d-th derivative, sum_k w_k f(x + k h) / h^d, has Taylor moments m_n = sum_k w_k k^n / n!. It is consistent when m_n = 0 for n < d and m_d = 1, and its order of accuracy is p when the first nonzero moment after m_d is m_{d+p} (LeVeque, Finite Difference Methods for ODEs and PDEs, sec. 1.5; Fornberg 1988). In the 2026-10 corpus sessions four repositories had FD tables to check this way, and one "8th-order" table had order 6.

stencil_moments(offsets, weights, upto)

Taylor moments m_0..m_upto of a stencil; exact when offsets and weights are integers or Fractions.

stencil_order(offsets, weights, deriv, *, tol=1e-10, max_order=20)

(scale, order): m_deriv (1 for a consistent stencil) and the order of accuracy, or (m_deriv, None) when a lower moment is nonzero (the stencil does not approximate this derivative).

check_stencil(offsets, weights, deriv, *, claimed=None, divisor=1, tol=1e-10)

A stencil approximates the deriv-th derivative, with the right scale, to the claimed order.

offsets are grid offsets such as [-1, 0, 1]; weights as written in the code, before division by divisor (pass 12 for weights written as integers over 12). Integer and Fraction weights are checked exactly. FAIL for a lower moment that is not zero, a scale other than 1 (a missing factor such as 1/2), or an order below claimed.

jacobi_moment(k, alpha=0.0, beta=0.0)

Exact integral of x^k (1-x)^alpha (1+x)^beta over [-1, 1], alpha, beta > -1.

m_0 = 2^(alpha+beta+1) B(alpha+1, beta+1); integrating d/dx [x^k (1-x)^(alpha+1) (1+x)^(beta+1)] = 0 gives (k + alpha + beta + 2) m_(k+1) = k m_(k-1) + (beta - alpha) m_k, free of the cancellation of the binomial expansion.

check_quadrature(nodes, weights, *, claimed, interval=(-1.0, 1.0), alpha=0.0, beta=0.0, moments=None, rtol=1e-10)

A quadrature rule integrates every polynomial up to degree claimed exactly.

An n-point Gauss rule is exact to degree 2n-1, Gauss-Radau to 2n-2, Gauss-Lobatto to 2n-3 (Abramowitz and Stegun 25.4; Hildebrand, Introduction to Numerical Analysis, ch. 8). The weight function is the Jacobi weight (1-x)^alpha (1+x)^beta on [-1, 1]; nodes on another interval are mapped there (alpha = beta = 0 only). For another weight pass moments(k), the exact integral of x^k on interval. Monomials are taken about the centre of the interval to keep the sums well conditioned. In the 2026-11 corpus session this found Radau weights that were wrong for alpha != beta.

physverify.fluxes

Exactness properties every numerical flux of a conservation law must have (Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., ch. 10-11; LeVeque, Finite Volume Methods, sec. 12.1):

• consistency F(W, W) = f(W), the physical flux, for any state W
• upwinding when every wave speed is positive on both sides, F(WL, WR) = f(WL) (all negative: f(WR)); Roe, HLL, HLLC and Rusanov fluxes with exact wave-speed bounds satisfy it
• mirror reflecting the normal velocity and swapping the sides gives F(R WR, R WL) = -R F(WL, WR), with R the reflection acting on the flux components

Each needs only the flux function and the physical flux written from the equations, not a solver run. In the 2026-10 corpus sessions these checks found both bugs in one Euler solver: a Roe flux that broke upwinding when density and tangential velocity both jumped, and an HLLE flux built for the wrong fluid.

check_flux(flux, physical, states, *, upwind_pairs=(), downwind_pairs=(), mirror=None, rtol=1e-10, atol=1e-12)

Consistency, upwinding and (optionally) mirror symmetry of a numerical flux flux(WL, WR).

physical(W) is the exact flux written from the equations; states are the states W for the consistency check. upwind_pairs are (WL, WR) with all wave speeds positive on both sides (F must equal f(WL)) and downwind_pairs all negative (F must equal f(WR)); the caller builds them from the sound speed. mirror = (reflect, signs): reflect(W) flips the normal velocity of a state and signs holds +1 or -1 per flux component (for the 1D Euler equations (-1, +1, -1): mass and energy fluxes change sign, the momentum flux does not); the check is F(reflect(WR), reflect(WL)) = signs * F(WL, WR), over all pairs.

physverify.operators

Invariants of discretised operators and solvers that need no exact solution.

Four properties a correct discretisation has whatever its order:

• check_nullspace: a discrete Laplacian, gradient or divergence applied to a constant field gives zero at every interior point (a constant is in the null space of the continuous operator, and consistent stencils keep it there).
• check_uniform_state: a uniform state with no flow and no source is a steady solution of conservation laws, so a time step must leave it unchanged.
• check_advection_direction: a feature carried by a velocity u moves by u t (Toro, Riemann Solvers, ch. 2). Magnitude checks pass for a solver that integrates U_t - F_x = 0 instead of U_t + F_x = 0; this one does not.
• check_force_gradient: a conservative force is minus the gradient of its potential, F = -dU/dr.

In the 2026-11-01 corpus sessions these four found four bugs that exact-solution comparisons had missed or needed far more work for: a scrambled polar-grid diffusion matrix (row sums 22 and 44 instead of 0), a Lax-Wendroff scheme that heated still water by 0.9 K in 5 s, an Euler solver that moved fluid at -u, and a screened-Coulomb force at 36-59 % of -dU/dr. Each passed on the correct codes it was tried on.

check_path tests the geometry a builder produced: a wire or path has no jumps and closes where it should (corpus session of 2026-11-06: a wire helper that always continued from vertex 1 gave a square coil 39 % too little field).

check_nullspace(operator, shape, *, field=None, interior=None, rtol=1e-10, scale=None)

An operator applied to a field in its null space (a constant by default) is zero away from the boundary.

operator is a matrix (dense or scipy sparse, acting on the flattened field) or a function field -> field. shape is the grid shape; field defaults to ones. interior is a tuple of slices, an int width, or None for one cell in from every edge (boundary rows hold boundary conditions, which need not vanish). The tolerance is relative to scale, by default the size of the operator applied to a non-constant field (so 1/dx^2 entries do not make round-off look like a failure).

check_uniform_state(step, state, *, n_steps=10, rtol=1e-10, atol=0.0)

A uniform state with no flow and no source is unchanged by time stepping.

step(state) -> new_state advances one step; state is the uniform initial state (any array, or a tuple of arrays). FAIL names the largest change. A scheme that returns an updated value where a flux or rate is expected (or the reverse), a source that is never zero, or a wrong sign in a boundary term shows up here.

check_advection_direction(evolve, x, u, t, *, initial=None, rtol=0.1)

A feature carried by velocity u moves by u t, in the direction of u.

evolve(initial, t) -> field runs the code under test on a periodic 1-D grid x (uniform) for time t and returns the advected field (density, tracer, a velocity component...). initial defaults to a smooth bump at one third of the domain. The shift is measured from the periodic cross-correlation, so u t should be well below half the domain. FAIL when the feature moves the wrong way (a sign error in the flux divergence or in a neighbour-index convention) or by more than rtol off u t.

check_force_gradient(potential, force, r, *, rel_step=1e-06, rtol=1e-05, atol=0.0)

A conservative force equals minus the derivative of its potential: F(r) = -dU/dr.

potential(r) and force(r) take a scalar or array of separations (wrap a function that returns (U, F) as lambda r: f(r)[0] and lambda r: f(r)[1]). The derivative is a central difference with step rel_step * r. Radial convention: positive F is repulsive. FAIL lists the separations where F differs.

check_path(vertices, *, segment_lengths=None, closed=False, rtol=1e-08)

A path built by a helper has the segments it was asked for: no jumps, and closed if it should be.

vertices has shape (n, dim). segment_lengths is one length or one per segment. A correct solver gives a wrong answer on a wrong path, so test the builder separately from the field or transport calculation.

stretched_grid(n, a=0.0, b=1.0, *, ratio=1.1)

n + 1 nodes from a to b whose spacings grow by ratio from one cell to the next.

Run the operator invariants (check_nullspace, check_uniform_state) on such a grid as well as a uniform one: a west-face coefficient built from the east spacing gives the right answer whenever the spacings are equal, so it is invisible on a uniform grid (corpus study: row sums 24 to 184 instead of 0 on a stretched grid, exactly 0 on a uniform one). For a second-derivative operator pass a linear field too: it is in the null space on any grid, and a wrong spacing shows up there even when a constant still gives zero.

Runge-Kutta tables and steppers

physverify.butcher

Order of a Runge-Kutta method from its Butcher table, by the order conditions.

A Runge-Kutta method (A, b, c) has order p exactly when b . Phi(t) = 1 / gamma(t) for every rooted tree t with at most p vertices (Butcher). This is algebra on the coefficients: it needs no test problem, cannot be fooled by a lucky cancellation or a superconvergent test equation, and works the same for explicit and implicit methods.

Lessons from testing other people's code: one package's order labels were wrong for 3 of its tables (one claimed 4 where the table satisfies only the order-3 conditions), and another's embedded weights for Tsitouras 5(4) did not even sum to 1. Order studies on test equations found both, but were slower to adjudicate than this check.

trees(order) cached

All rooted trees with order vertices, each a sorted tuple of child subtrees.

order_residuals(A, b, max_order=6)

{order: largest |b . Phi(t) - 1/gamma(t)| over the trees of that order}.

tableau_order(A, b, *, max_order=6, tol=1e-08)

Highest p <= max_order such that every condition up to order p holds within tol.

linear_residuals(A, b, max_order=8)

{k: |b . A^(k-1) 1 - 1/k!|}: how far the stability function R(z) is from e^z at each power of z.

linear_order(A, b, *, max_order=8, tol=1e-08)

Order on linear constant-coefficient problems (y' = λy): highest p with R(z) = e^z + O(z^(p+1)).

It is at least the classical order and can be higher; a convergence test on y' = λy measures this number, not the classical one (Hairer & Wanner, Solving ODEs II, section IV.3).

stage_order(A, c=None, b=None, *, max_order=8, tol=1e-08)

Highest q with A c^(k-1) = c^k / k for k = 1..q (simplifying assumption C(q)).

With b given, b . c^(k-1) = 1/k (assumption B(q)) is required as well, the usual definition of stage order. A low stage order is what causes order reduction on stiff problems.

stability_polynomials(A, b)

(P, Q) coefficient arrays, lowest power first, with R(z) = P(z) / Q(z).

P(z) = det(I - zA + z 1 b^T) and Q(z) = det(I - zA) have degree at most s; they are recovered from their values on the unit circle, which is well conditioned.

stability_function(A, b)

R(z) = 1 + z b^T (I - zA)^(-1) 1, as a function of a complex z.

r_infinity(A, b)

lim R(z) as |z| -> infinity: 0 for an L-stable method, inf when R is a polynomial (explicit).

a_stability_margin(A, b)

min over real y of |Q(iy)|^2 - |P(iy)|^2, scaled by the largest coefficient; >= 0 on the imaginary axis means |R(iy)| <= 1. Returns -inf when R has a pole in the closed left half-plane or grows at infinity.

A method is A-stable exactly when the margin is >= 0 (up to round-off) and every pole of R lies in Re z > 0 (Hairer & Wanner II, section IV.3, the E-polynomial).

is_stiffly_accurate(A, b, *, tol=1e-12)

b equals the last row of A, so the step result is the last stage (Hairer & Wanner II, IV.6).

square_tableau(A, b, c=None)

(A, b, c) with A made square: a table stored as s x (s - 1) (the last column of an explicit table is zero and is often left out) or with an extra FSAL stage that has weight 0 in b is padded with zeros. Explicit FSAL tables (Dormand-Prince, Bogacki-Shampine) then check like any other.

symplectic_residual(A, b)

max |b_i a_ij + b_j a_ji - b_i b_j|: zero exactly for a symplectic Runge-Kutta method (Sanz-Serna 1988; Hairer, Lubich and Wanner, Geometric Numerical Integration, Thm. VI.4.3). Gauss methods are symplectic; no explicit Runge-Kutta method is.

check_tableau(A, b, c=None, *, claimed, label='', tol=1e-08, max_order=None, claimed_linear=None, claimed_stage=None, a_stable=None, l_stable=None, stiffly_accurate=None, embedded=None, claimed_embedded=None, symplectic=None)

PASS when the table has at least the claimed order and every other stated property; FAIL otherwise.

The default tolerance 1e-8 accepts coefficients typed with 9 or 10 digits (a real package had residuals of 2e-10 for that reason), while a missing order shows residuals of 1e-4 or more. The message gives the residuals, so a table whose conditions hold only to about 1e-10 is visible, and a claim that is lower than the true order is noted rather than failed. If c is given, the row sums of A must equal it (the order conditions assume they do, so a mismatch fails).

Optional claims, each checked only when given: claimed_linear order on y' = λy (the order a Dahlquist convergence test measures); see linear_order claimed_stage stage order, C(q) and B(q); see stage_order a_stable |R(z)| <= 1 on Re z <= 0; see a_stability_margin l_stable A-stable and R(inf) = 0 stiffly_accurate b equals the last row of A symplectic b_i a_ij + b_j a_ji = b_i b_j; see symplectic_residual True means "the method has it" and False means "the method does not"; a mismatch either way fails. embedded the weights of the embedded (error-estimate) method; with claimed_embedded, its order is checked too, and the error weights b - embedded must sum to 0. A table stored as s x (s - 1), or with an FSAL stage of weight 0, is padded first (square_tableau).

compare_tableau(A, b, name, *, tol=1e-10)

A code's copy of a published table (one of TABLES) has the published coefficients.

physverify.steppers

Reference time steppers and a fixed-step driver.

These are used (a) to build independent reference solutions and (b) as the stepper when only a right-hand side, not a whole solver, is under test.

integrate_fixed(f, y0, t0, t1, n_steps, stepper=rk4)

Advance y' = f(t, y) from t0 to t1 in n_steps equal steps.

References and constants

physverify.references

Exact reference solutions that testing public codes needed again and again.

Each one was written by hand in several adapters before it moved here (repository IDs as in results/corpus/SUMMARY.md). All are numpy-only and vectorised where that is natural.

ising2d_exact(L, T, J=1.0, h=0.0)

Exact thermal averages of the periodic L x L Ising model by enumerating all 2^(L*L) states.

H = -J sum_ s_i s_j - h sum_i s_i, k_B = 1. Returns per-site , <|m|>, , , the specific heat per site and Z. Practical up to L = 4 (65536 states); used for D04, H06, H55.

langevin(x)

L(x) = coth x - 1/x, with the series near 0.

heisenberg_open_chain(N, J, T)

Classical unit-vector Heisenberg chain with open ends, H = -J sum S_i.S_{i+1}: exact per-site energy and specific heat (each bond is independent). Used for H64.

heat_sine_mode(x, t, L=1.0, kappa=1.0, mode=1)

u_t = kappa u_xx on [0, L], u = 0 at both ends, u(x, 0) = sin(mode pi x / L).

advection_diffusion_gaussian(x, t, x0, sigma0, v=0.0, D=0.0)

Gaussian of peak 1 and width sigma0 at x0, carried at speed v and spread by D on the whole line: sigma^2 = sigma0^2 + 2 D t, peak sigma0/sigma. Used for H25.

semi_infinite_step(x, t, alpha, T0, Ts)

Semi-infinite solid at T0, surface held at Ts from t = 0: T = Ts + (T0 - Ts) erf(x / 2 sqrt(alpha t)).

double_pendulum_acc(th1, w1, th2, w2, m1, m2, l1, l2, g=9.81)

Angular accelerations of a double pendulum (point masses, massless rods, angles from the downward vertical), from the Lagrangian mass matrix. Used for D01, E01, H03, H27.

cantilever_tip(P, L, E, I)

Euler-Bernoulli cantilever with an end load: tip deflection P L^3/(3EI) and rotation P L^2/(2EI).

infinite_well_energy(n, W, hbar=1.0, m=1.0)

E_n = n^2 pi^2 hbar^2 / (2 m W^2), n = 1, 2, ... Used for H19.

barrier_transmission(E, V, a, hbar=1.0, m=1.0)

Transmission probability through a square barrier of height V and width a (E != V). Used for H40.

tfim_open_ground_energy(N, J, h)

Ground energy of the open transverse-field Ising chain H = J sum Z_i Z_{i+1} + h sum X_i, from the free-fermion (Jordan-Wigner, Bogoliubov-de Gennes) solution; the sign of J does not matter. Used for H22.

morse_levels(v, omega_e, omega_e_x_e)

Morse vibrational levels above the potential minimum, G(v) = w_e (v + 1/2) - w_e x_e (v + 1/2)^2, exact for the potential D (1 - exp(-a (r - r_e)))^2 with D = w_e^2 / (4 w_e x_e). Used for J222, J414.

morse_turning_width(E, D, a)

r_max - r_min of the Morse potential D (1 - exp(-a (r - r_e)))^2 at energy E (0 < E < D): (1/a) ln((1 + q)/(1 - q)) with q = sqrt(E/D). An RKR inversion reproduces this exactly from G(v) alone. Used for J414.

schwarzschild_qnm(l, s=2)

Fundamental (n = 0) quasinormal frequency M omega of a Schwarzschild black hole for a field of spin s (0 scalar, 2 gravitational) and multipole l; 5-digit literature values (Leaver 1985; Berti, Cardoso and Starinets 2009). A ringdown that matches s = 0 when the code claims gravitational waves points to the potential term (1 - s^2) 2M/r^3. Used for J216.

kerr_isco(a, prograde=True)

Radius of the innermost stable circular orbit of a Kerr black hole (Bardeen, Press and Teukolsky 1972), in units of M; a is the dimensionless spin, 0 <= a <= 1. Gives 6 at a = 0, and 1 or 9 at a = 1. Used for J116.

schwarzschild_deflection(b, nodes=200)

Exact deflection angle (radians) of light with impact parameter b > 3 sqrt(3) M past a Schwarzschild black hole, G = c = M = 1: 2 int_0^{u0} du / sqrt(1/b^2 - u^2 + 2 u^3) - pi, u0 the turning point. The substitution u = u0 (1 - s^2) removes the endpoint singularity; Gauss-Legendre on s. Tends to 4/b. Used for J243.

inhour_stable_period(rho, Lambda, beta_i, lambda_i)

Asymptotic (stable) reactor period in seconds from the inhour equation rho = w Lambda + sum_i w beta_i / (w + lambda_i): the largest root, positive for rho > 0 and in (-min lambda_i, 0) for rho < 0 (it tends to -1/min lambda_i, about -80 s for U-235, after a large negative step). Used for J192.

zernike_osa(n, m, rho, theta)

Zernike polynomial Z_n^m in the OSA/ANSI convention (ANSI Z80.28; Thibos et al. 2002): N R_n^|m|(rho) times cos(m theta) for m >= 0 and sin(|m| theta) for m < 0, with N = sqrt(2 (n + 1) / (1 + delta_m0)), so that the mean of Z^2 over the unit disk is 1. Noll's ordering uses the same normalisation. Two corpus repos needed it (L173 right, L208 normalised as sqrt(n + 1) for every m).

physverify.constants

Find numeric literals in source code that are close to, but not equal to, a physical constant.

Lessons from testing other people's code: a gravitational constant in au^3 / (Msun yr^2) of 39.4769 (Gaussian, 365.25-day year) next to an orbital speed of 2 pi (which assumes 4 pi^2 = 39.4784), g = 9 in a projectile game, pi typed as 3.14. None of these is necessarily a bug, but each is worth one look, and an adapter's reference must use the code's value, not the textbook one.

scan(path) returns notes, never failures: the code may be right to round. Several CODATA editions are accepted. When writing a reference, use the edition the code says it uses: one adapter failed at 1.35e-9 by taking the 2022 proton mass where the code stated 2018.

check_value(value, *, max_rel=0.1)

(constant, nearest accepted value, relative difference) for a near miss, else None.

A value equal to an accepted value, or equal to it to the precision written with at least 4 digits, is not a near miss: 3.1416 for pi and 9.81 for g are fine; 3.14 and 9 are noted.

scan(path, *, max_rel=0.1)

Near-miss constants among the numeric literals of a Python file.

Literals with 3 or more significant digits are compared with every constant. Integers and shorter literals (9, 3.1, 10) are so common that they are noted only when assigned to a name that suggests the constant (g = 9).

Statistics and Monte Carlo

physverify.stochastic.stats

Statistics for Monte Carlo checks: error bars, autocorrelation and verdicts.

Design (docs/plan_ising.md, section 4):

• Every estimate is a mean over independent chains (each chain is burned in and then averaged). Chains are independent by construction, so the standard error is the spread of the chain means over sqrt(n_chains), with no autocorrelation model needed. Within one chain, tau_int and blocked_error give the error bar for codes that report one.
• A check makes m comparisons with exact values. Each uses a two-sided Student-t test with n_chains - 1 degrees of freedom at level alpha / m (Bonferroni), so the false-alarm rate of a whole check is at most alpha.
• Verdicts: FAIL when some comparison is significant; INCONCLUSIVE when none is, but the test could not have detected a deviation as large as the stated resolution (critical value times standard error > resolution); PASS otherwise.

t_sf(x, nu)

P(T > x) for Student's t with nu degrees of freedom (x >= 0), by Simpson's rule on the pdf.

t_critical(alpha_two_sided, nu) cached

x with P(|T| > x) = alpha, by bisection (cached: checks ask for the same few values).

tau_int(series, c=6.0)

Integrated autocorrelation time with Sokal's automatic window (window W >= c * tau).

blocked_error(series)

Standard error of the mean of a correlated series: sqrt(2 tau var / n).

compare_chain_means(name, per_chain_values, exact, resolution, min_se=0.0)

Mean over independent chains, with the spread of the chain means as its standard error.

min_se is a floor for the standard error: for the frequency of a rare event the chain-to-chain spread can be zero when the event never happens, which would make any small expected probability look like a certain deviation.

verdict(comparisons, alpha=ALPHA_CHECK)

Bonferroni-corrected t-tests; see the module docstring for PASS, FAIL and INCONCLUSIVE.

detectable_deviation(atol, std, n, *, power=0.9)

Smallest true deviation of a mean that the test |mean - exact| <= atol fails with probability power.

The sample mean has standard error std / sqrt(n), so a deviation d is caught with probability power when d >= atol + z_power * std / sqrt(n) (normal approximation). A moment test whose tolerance comes from the sample size can pass while a term is wrong: in one stochastic-integrator test atol was 0.078 and a deviation of 0.049 from a wrong term passed. Report this number next to such a test.

check_moment(values, exact, atol, *, resolution=None, power=0.9)

The mean of independent values is within atol of exact, with the size of deviation the test can see.

FAIL when |mean - exact| > atol. INCONCLUSIVE when it passes but a deviation of resolution (the smallest error that matters for the physics) would pass too with probability above 1 - power. PASS otherwise.

physverify.stochastic.checks

Checks for an Ising sampler (see physverify.stochastic.sampler for the interface).

Deterministic (plain PASS / FAIL, no statistics):

bookkeeping tracked E and M, the code's own measurement, and an independent recomputation from the lattice all agree exactly after every sweep scripted_rng the code's sweeps, fed recorded random numbers, give exactly the lattices of the trusted reference replaying the same numbers reproducibility the same seed gives the same run; different seeds do not

Statistical (PASS / FAIL / INCONCLUSIVE, false-alarm rate per check stated in the message):

exact_distribution fractions of time at each energy level and each |m| level, L = 4, against exact enumeration moments , , <|m|>, at temperatures below, at and above T_c (L = 4, enumeration), in a field, and , at L = 8 (Kaufman's closed form) acceptance single updates at controlled sites: acceptance frequency for each energy change against min(1, exp(-dE / T)) (Metropolis) or the heat-bath probability symmetry translation and spin-flip symmetry: every bond's mean correlation and every site's mean spin against the exact values equilibration hot and cold starts agree from the first reported sweep error_bars the code's reported error bars against the actual spread of independent chains

Each check uses its own fixed seed, so the whole suite is deterministic.

RecordingRNG

A numpy Generator that records every integers() and random() call, in order.

applicable(sampler, family)

False for a check the sampler cannot support: missing methods, or not a single-site algorithm.

ising_checks(sampler, **overrides)

The checks that apply to this sampler (see applicable).

overrides maps a check name to keyword arguments for it, for example {"moments": {"n_chains": 128}} for a slow sampler. Called directly, a check that does not apply reports INCONCLUSIVE instead.

physverify.stochastic.directions

Direction sampling in Monte Carlo transport codes (issue #47).

A sampler returns unit vectors; mu = d . n is the cosine to a normal n. The exact laws:

• isotropic on the sphere: mu uniform on [-1, 1], CDF (mu + 1) / 2;
• isotropic on the hemisphere about n: mu uniform on [0, 1], CDF mu;
• Lambertian (cosine-weighted, diffuse wall emission or reflection): density 2 mu on [0, 1], CDF mu^2, mean 2/3 (Modest, Radiative Heat Transfer, sec. 1.6). Sampling mu = sqrt(xi) gives it; mu = xi gives the hemisphere law.

The azimuth about n must be uniform for all three. The test is a Kolmogorov-Smirnov test on mu and one on the azimuth, each at alpha / 2, so the false-alarm rate of the check is at most alpha (measured on correct samplers at alpha = 1e-3, n = 20 000: 0 of 2000 runs for the Lambert law, 1 of 2000 for the isotropic law). In the corpus study a phonon-transport code drew wall reflections with mean mu 0.700 instead of 2/3; at the default n = 200 000 that deviation is found.

ks_pvalue(d, n)

Asymptotic p-value of the one-sample Kolmogorov-Smirnov statistic d for n samples (Stephens' small-sample correction; Numerical Recipes 3rd ed., sec. 14.3.3).

check_direction_sampler(sample, *, law='isotropic', normal=(0.0, 0.0, 1.0), n=200000, seed=0, alpha=0.001)

Directions from sample(rng, n) -> (n, 3) follow law ("isotropic", "hemisphere" or "lambert") about normal.

FAIL when a vector is not of unit length, when mu = d . normal does not follow the law, or when the azimuth about the normal is not uniform. The message gives the sample mean of mu next to the exact one.

physverify.stochastic.ito

Exact properties of iterated Ito integrals, for codes that generate them for Milstein or order 1.5 schemes.

From Ito calculus (Kloeden and Platen, Numerical Solution of SDEs, ch. 5), on a step h with W = W(h): I_ii = (W_i^2 - h)/2; I_ij + I_ji = W_i W_j for i != j; I_iii = (W_i^3 - 3 h W_i)/6; and by the Ito isometry E[I_ij^2] = h^2/2, E[I_ijk^2] = h^3/6 for every index choice.

Generators with several noises use a truncated series, so an identity may hold only as the series grows: a truncation error falls with the series length and a wrong formula does not. In the corpus session of 2026-11-06 one generator had an I_111 residual of 0.079 and 0.081 RMS at series lengths 10 and 20, and E[I_ijk^2] from 0.116 to 0.206 instead of 1/6. False alarms on correct integrals are measured in tests/test_new_checks_2026_11.py and stated in docs/checks/numerical-controls.md.

check_ito_integrals(generate, ps=(None,), *, h=1.0, rtol=0.02, margin=0.05, z=5.0)

generate(p) returns (W, I2, I3): samples of shape (n, m), (n, m, m), (n, m, m, m); I2 or I3 may be None.

ps are series lengths in increasing order. With one entry, every identity must hold within rtol h^k RMS. With two or more, an identity above rtol at the longest series must have fallen by at least 25% from the shortest (the error of a truncated Levy-area series goes as 1/sqrt(p) or faster, so use lengths a factor of 4 or more apart). The isometry at the longest series allows margin relative plus z standard errors.

Loading code

physverify.loaders

Import a module from a path without touching its folder.

Used by the adapters that test other people's code. Nothing is written next to the module (no pycache), and optional plotting imports can be stubbed so a simulation file can be imported on a machine without matplotlib.

Placeholder

Stand-in for a name from a module that is not installed.

Calling it returns another placeholder, except that calling it with a single function returns that function, so stubbed decorators (@app.route(...), @jit) leave the decorated code usable. Attribute access also returns a placeholder. A tested code path that really uses the stub gets a placeholder back instead of a result, which the checks then report.

stub_modules(spec, *, force=False)

Install stub modules for imports a tested code path does not use.

spec is an iterable of module names, or a mapping from a module name to the names it must provide ({"sklearn.metrics": ["accuracy_score"]}), so that from sklearn.metrics import accuracy_score succeeds. Parent packages are created as needed. A module that imports for real is left alone unless force is true. Returns the names installed.

run_script(path, *, argv=(), stubs=None, inputs=None)

Run a whole script as main and return (its globals, what it printed).

For scripts whose results are only module-level variables (coursework, notebooks exported to .py, demos). Matplotlib is switched to a non-interactive backend and plt.show() does nothing, so the run never waits for a window.

inputs answers the script's input() prompts in order (interactive scripts); a prompt with no answer left raises EOFError instead of waiting.

load_definitions_from_path(name, path, *, keep=(), inner=(), namespace=None, stubs=())

Load only the definitions of a script whose top level runs a simulation.

Many research scripts parse argv, open a window or run a long simulation at import. This executes only the top-level imports, functions and classes, and the top-level assignments to the names in keep (constants the functions use). inner names functions defined inside other functions (for example marimo or Jupyter-exported cells); they are lifted to module level, and their imports must then come from namespace.

Nothing else in the file runs, so a function that relies on a global set by the skipped code fails with NameError; set it on the returned module.

stubs names modules (or maps module names to the names they must provide, see stub_modules) to replace with stubs when they are not installed (pygame, IPython, orbax, ...), for imports the tested functions do not use.

extract_functions(path, names, *, namespace=None, stubs=())

Load the named top-level functions or classes of one file, plus the same-file definitions they use.

For a package that cannot be installed or imported as a whole (a compiled extension, a dependency that needs a newer NumPy, an import of a GUI toolkit), when the function under test is plain Python. The named definitions are taken with every top-level function, class or simple assignment they refer to, followed transitively (a curvature that calls centered_diff brings it along). Each top-level import is run on its own; one that fails is skipped and listed in the returned module's __skipped_imports__, so a NameError later points at it. Names the code needs from other files of the package can be passed in namespace (for example a stub gf object with the one function used).

Nothing else in the file runs. The returned module records what was loaded in __extracted__.

load_repaired(name, path, fix, *, stubs=())

Import a file that does not run as published, after a recorded repair of its source.

fix is a dict {old text: new text} (each old must occur in the file, else ValueError) or a function source -> source. The repaired source is compiled in memory with the original file name, so tracebacks point at the right lines; nothing is written to disk. The module's __repairs__ lists what was changed, and repairs_note(module) gives a sentence to put in the check's message, so that a PASS on repaired code is never reported as a PASS on the original (batch 7 had scripts that did not parse, one with a stray line, one written for Python 2).

repairs_note(module)

'tested after repairs: ...' for a module from load_repaired, or '' when nothing was changed.