API reference

Fitting

pysnape.fit(model, data, *, unknowns=None, known=None, constraints=None, knots='auto', order=None, multistart=True, weights=None, solver='auto', verbose=False)

Estimate the unknown parameters of an ODE or PDE from data.

Parameters:

Name	Type	Description	Default
model	str | Sequence[str] | Model	The equation as written on paper ("x_tt + c*x_t + k*x = 0"), a list of equations sharing parameters, or a Model such as ps.models.burgers(). Fields, coordinates and known inputs are read from data.	required
data	Data	The measurements, as Data.	required
unknowns	Mapping[str, float] | Sequence[str] | None	Starting values {"c": 1.0, ...}, or just the names (each starts at 1). Required with equation strings; with a Model it overrides the model's starting values.	None
known	Mapping[str, float] | None	Fixed constants and sign constraints (">=0", "<=0"), as in build_model; only with equation strings.	None
constraints	Mapping[str, float] | None	Fixed constants and sign constraints (">=0", "<=0"), as in build_model; only with equation strings.	None
knots	Knots | Sequence[Knots]	"auto" (about one knot per three grid points), one count for all coordinates, a mapping per coordinate, or a list of such candidates: each is fitted and the one with the lowest data misfit is kept.	'auto'
order	Orders	B-spline order per coordinate; by default the highest derivative + 3.	None
multistart	bool	Also search from a grid of data-scaled starting points at every weight (slower, safer on noisy data). Turn off for clean data or when the start is good.	True
weights	Sequence[float] | None	Equation weights for Nested, in balanced units; the default is Nested's. Fewer, smaller weights are faster.	None
solver	str	Sparse solver for the spline fit ("auto" uses CHOLMOD when installed).	'auto'
verbose	bool	Print each basis tried and each weight stage as it finishes, with the elapsed time: useful on large grids, where a fit can take many minutes.	False

Returns:

Type	Description
Fit

pysnape.Fit dataclass

What fit found.

theta holds the estimates; field() the fitted fields and their derivatives on the data grid; summary() a readable table; plot() measured against fitted.

data_misfit property

||H beta - y|| / ||y|| at the observed points.

equation_residual property

||g|| / ||g_0||: how far the fitted fields are from satisfying the equation.

field(name=None, deriv=None)

A fitted field, or its derivative (deriv="xx"), as an array on the data grid.

fitted(name=None)

The fitted counterpart of what the data holds for a field (its observed derivative).

residual()

The equation residual on the grid; shape (n_equations, *grid) when several.

refit(data, *, start=None, multistart=None)

Fit new data with this fit's model, basis and settings.

start defaults to this fit's starting values; multistart to its setting.

summary()

A plain-text report: the equations, the basis, the estimates and how well it fits.

plot(name=None)

Measured against fitted for one field; returns the matplotlib axes.

1-D data: both curves, on one Axes. 2-D data: measured, fitted and their difference as images, on an array of three. axes.figure (or axes[0].figure) is the figure. Returning the axes, as pandas and seaborn do, keeps a notebook from drawing the figure twice. Needs matplotlib (pip install "pysnape[plot]").

pysnape.fitting.basis_spec(model, data, knots='auto', order=None)

The {coord: {"knots": k, "order": o}} spec fit uses.

knots is "auto", one count for every coordinate, or a mapping per coordinate (missing coordinates get "auto"). order is None (auto), one order for every coordinate, or a mapping (missing coordinates get the auto order).

Uncertainty

pysnape.ensemble(fit, n=20, *, noise=None, snr_db=None, seed=0, jobs=1)

Refit n noisy copies of the fit's (noise-free) data.

Draw i adds noise of noise percent (or snr_db) with seed seed + i, the papers' convention, and refits from the fit's starting values with its settings.

pysnape.bootstrap(fit, n=100, *, seed=0, multistart=False, jobs=1)

Residual bootstrap: refit the fitted fields plus resampled residuals, n times.

Residuals are resampled with replacement within each field, at the observed points only. Each refit starts from the fit's estimate, without multi-start unless asked.

pysnape.Ensemble dataclass

Parameter estimates over repeated fits.

thetas is (n, n_params); mean, std and cov (coefficient of variation, in percent) summarise it, and interval() gives percentile bounds.

cov property

Coefficient of variation, 100 * std / |mean|, in percent.

interval(level=0.95)

Percentile interval holding level of the estimates, per parameter.

summary(truth=None, level=0.95)

A table of mean, COV and interval; with truth, the bias of the mean too.

Data

pysnape.Data dataclass

Fields, coordinates and known inputs on a grid.

Parameters:

Name	Type	Description	Default
fields	Mapping[str, ArrayLike]	{"u": array} -- what was measured of each field, shaped like the grid.	required
coords	Mapping[str, ArrayLike]	{"x": x, "t": t} -- 1-D, strictly increasing; the key order is the coordinate order and fixes the array layout.	required
inputs	Mapping[str, ArrayLike]	Known signals on the grid, by name.	dict()
observe	Mapping[str, str]	{"x1": "tt"} -- the derivative suffix that names what the array under a field holds; absent means the field itself.	dict()
mask	ArrayLike | None	Boolean array shaped like the grid: the grid points that were observed. None means all of them.	None

observed_rows property

Indices (in vector order) of the observed grid points; None if all.

vector(name)

A field or input as a vector, first coordinate fastest (MATLAB reshape(u, [], 1)).

array(vec)

Inverse of vector.

grid()

Each coordinate broadcast over the grid (np.meshgrid(..., indexing="ij")).

replace(**changes)

A copy with some attributes replaced (re-validated).

add_noise(snr_db=None, *, percent=None, seed=None, fields=None)

A copy with white Gaussian noise on the fields, awgn(..., 'measured') style.

One generator, seeded once, serves every field in turn, so seed=100 names one noise realisation of the whole dataset -- the thing an ensemble loops over. Inputs are known signals and are left alone.

subsample(fraction, *, seed=None, per_axis=True)

A copy observing a random fraction of the grid (the compressive-sensing setup).

per_axis=True keeps a random sub-grid -- fraction of the x samples times fraction of the t samples, the papers' "70 % x 70 %"; False keeps random points. The equation is still enforced everywhere; only the data term shrinks.

save(path)

Write everything to one .npz; load reads it back.

load(path) classmethod

Read a file written by save.

from_file(path, *, fields, coords, inputs=None, observe=None, mask=None) classmethod

Build from a file whose array names differ from ours: fields={"u": "usol"}.

pysnape.data.noise.add_noise(x, snr_db=None, *, percent=None, seed=None, rng=None)

awgn(x, snr_db, 'measured'): x plus white Gaussian noise at the given SNR.

The noise power is the measured power of x (mean of |x|**2) divided by 10**(snr_db/10). Pass seed for a fresh generator, or rng to keep drawing from an existing one -- the fields of one dataset share a generator that way.

pysnape.data.datasets

Published datasets, fetched on first use: load("burgers") gives a Data.

Each dataset's source, citation, license note, checksum and layout live in manifest.yaml next to this module. Files are downloaded from their source into a cache directory -- $PYSNAPE_DATA if set, else ~/.cache/pysnape -- and checked against their sha256 before use; nothing is redistributed with pysnape. Cite the source paper when you use one (info(name)["citation"]).

available()

Names load accepts.

info(name)

Description, equation, source, citation, license note, url, checksum and size.

cache_dir(directory=None)

Where downloads go: directory, else $PYSNAPE_DATA, else ~/.cache/pysnape.

fetch(name, directory=None)

The local path of a dataset's file, downloading and verifying it if needed.

load(name, directory=None)

A dataset as Data, with field u on coordinates ("x", "t").

Fields stored as complex with negligible imaginary parts (spectral-solver round-off) are returned real. Pair it with the matching model: ps.models.<info(name)["model"]>().

Models

pysnape.build_model(equations, fields, coords, unknowns=None, inputs=(), known=None, constraints=None)

Parse one or more equations into a Model with shared unknowns.

A parameter name that appears in several equations is the same parameter -- that is what makes the MDOF chain (k2 in equations 1 and 2) work with no special casing.

pysnape.model.zoo

Model builders: the systems of the two papers as ready-made Model objects.

A builder does no numerics. It is build_model with the equation, the field and coordinate names, the starting point and the sign constraints of one paper case already filled in, so:

model = zoo.burgers()

is the same object as:

build_model(
    "u_t + theta1*u*u_x + theta2*u_xx = 0", fields=["u"], coords=["x", "t"],
    unknowns={"theta1": 3.0, "theta2": 3.0}, constraints={"theta1": ">=0"},
)

Conventions shared by every builder:

• The unknowns are called theta1, theta2, ... in the order the papers list them, so NestedResult.theta_array lines up with TRUTH.
• Equations carry the signs that give the parameters the values the papers report: KdV is u_t + 6 u u_x + u_xxx = 0, giving (6, 1).
• theta0 is a default starting point. constraints="paper" is the sign constraint the papers used for that case; None means none, and a mapping or a ">=0" / "<=0" string is passed to build_model unchanged.
• The data must use the same field, input and coordinate names in the same order, and the basis in each coordinate needs order model.max_deriv(c) + 2 (the papers' rule; the assembler enforces it). PDEs use the field u on ("x", "t"); the ODEs use x on ("t",).

burgers(theta0=None, constraints=PAPER)

u_t + theta1 u u_x + theta2 u_xx = 0, truth (1.0, -0.1).

Starts at (3, 3); keeps theta1 >= 0.

kdv(theta0=None, constraints=PAPER)

u_t + theta1 u u_x + theta2 u_xxx = 0, truth (6, 1).

Starts at (3, 3); keeps both parameters >= 0. The x basis needs order 5.

ks(theta0=None, constraints=PAPER)

Kuramoto-Sivashinsky, u_t + theta1 u u_x + theta2 u_xx + theta3 u_xxxx = 0.

Truth (1, 1, 1); starts at (-1, -1, -1); no constraints. The x basis needs order 6.

duffing(theta0=None, constraints=PAPER)

Forced Duffing, x_tt + theta1 x_t + theta2 x + theta3 x^3 = p, truth (0.5, -1, 1).

p is the known forcing, 0.42 cos t in the papers, supplied as the data input p (pysnape.data.generators.ode.duffing returns it). Starts at (5, 5, 5); no projections.

vanderpol(theta0=None, constraints=PAPER)

x_tt + theta1 x_t + theta2 x^2 x_t + theta3 x = 0, truth (-mu, mu, 1) = (-8, 8, 1).

Starts at (5, 5, 5); no constraints. x^2 x_t freezes x twice and keeps x_t live.

wave2d(theta0=None, constraints=PAPER)

2-D wave, u_tt - theta1 u_xx - theta2 u_yy = 0 on (x, y, t), truth (1, 1).

The papers report (1, 1): the wave coefficient of u_tt = c (u_xx + u_yy). Starts at (5, 5). No projections. Data: pysnape.data.generators.pde2d.wave.

mdof_chain(n=5, masses=1.0, theta0=None, constraints='>=0')

A chain of n masses fixed at one end: unknown springs k1..kn, dampers c1..cn.

Equation i is the force balance on mass i:

m_i x_i'' + c_i (x_i' - x_{i-1}') - c_{i+1} (x_{i+1}' - x_i')
          + k_i (x_i - x_{i-1})   - k_{i+1} (x_{i+1} - x_i)   = f_i

with x_0 = 0 (the base) and no i + 1 terms on the last mass. The masses are known; f1..fn are known force inputs (zero where a mass is not driven). Every k_j and c_j appears in two equations, which the solver handles as one shared parameter. Starts at k = 5, c = 0.5 (theta0 overrides any of them); all parameters are kept >= 0. Data: pysnape.data.generators.ode.mdof_chain, displacements or accelerations.

Generators

pysnape.data.generators.ode

Synthetic ODE data: the papers' Duffing and Van der Pol series and an MDOF chain.

Each generator integrates the system with an 8th-order Runge-Kutta method (DOP853, rtol=1e-12, atol=1e-13), so that any error in a fit comes from the noise and the method, not from the integrator.

The default is DOP853 rather than LSODA on the strength of a measurement (scipy 1.18, 2026-09-20). DOP853 at 1e-12 and 1e-13 and Radau at 1e-12 agree with each other to 6e-10 over Duffing's first 60 time units and to 4e-6 over the full 200, and to 4e-11 on Van der Pol (mu = 8, 50 time units). LSODA at rtol=1e-10 is off from all three by 2.7e-6 at 60 and by 5e-3 at 200 on Duffing, comparable to the noise of the 0.5 % experiments, and by 7e-7 on Van der Pol. Its tolerance is not its global error.

Only the displacement x is returned as a field -- the velocity is the integrator's state, not something the papers observe. Noise is a separate step, data.add_noise.

The defaults are the papers' systems, so with the default arguments the true parameters are the ones in pysnape.model.zoo.TRUTH:

Duffing      x'' + delta x' + alpha x + beta x^3 = amp cos(omega t)
             theta = (delta, alpha, beta)
Van der Pol  x'' - mu (1 - x^2) x' + x = 0
             theta = (-mu, mu, 1)

mdof_chain is not from the papers: a chain of masses, springs and dampers fixed at one end, driven by known forces, observed as displacements or accelerations.

duffing(t_end=200.0, n=4001, *, delta=0.5, alpha=-1.0, beta=1.0, amp=0.42, omega=1.0, x0=(0.5, 0.1), method='DOP853', rtol=1e-12, atol=1e-13)

The forced Duffing oscillator of the SNAPE paper.

Solves x'' + delta x' + alpha x + beta x^3 = amp cos(omega t). The default grid is 0 to 200 in steps of 0.05 (4001 points). The forcing is returned as the known input p = amp cos(omega t), which is what zoo.duffing() expects on the right-hand side; the true parameters are (delta, alpha, beta).

vanderpol(t_end=50.0, n=5000, *, mu=8.0, x0=(2.0, 0.0), method='DOP853', rtol=1e-12, atol=1e-13)

The Van der Pol oscillator of the SNAPE paper.

Solves x'' = mu (1 - x^2) x' - x. The default grid is 5000 points from 0 to 50, with mu = 8: a relaxation oscillation with sharp fronts. The true parameters are (-mu, mu, 1).

chain_matrices(m, k, c)

(M, C, K) of a chain fixed at its first spring and free at its last mass.

multisine(t, *, low=0.2, high=8.0, count=40, amplitude=1.0, seed=0)

A band-limited force: count sines from low to high rad/s, random phases.

Scaled so its RMS is amplitude. Exciting the whole band keeps every mode of the structure visible in the response, which is what makes all the parameters identifiable.

mdof_chain(n_dof=5, t_end=40.0, n=2001, *, m=1.0, k=None, c=None, forced=None, amplitude=1.0, seed=0, observe='displacement', method='DOP853', rtol=1e-12, atol=1e-13)

A fixed-free chain of n_dof masses, driven by known forces, at rest at t = 0.

Spring k_i and damper c_i join mass i - 1 to mass i, mass 0 being the fixed base, so M x'' + C x' + K x = f with tridiagonal C and K (chain_matrices); pysnape.model.zoo.mdof_chain writes the same equations.

Defaults: unit masses, k falling from 12 to 4 and c from 0.6 to 0.2 along the chain (for five masses k = (12, 10, 8, 6, 4), c = (0.6, 0.5, 0.4, 0.3, 0.2)), and one multisine force on the free end -- a shaker on the top floor. forced lists the masses (1-based) that are driven, each with its own seed.

Returns fields x1 .. xn: displacements, or with observe="acceleration" the accelerations (computed from the equations of motion, not differenced), marked with Data.observe so that fit knows what was measured. Inputs f1 .. fn hold the forces, zero on the masses that are not driven.

pysnape.data.generators.pde2d

Synthetic 2-D PDE data: the 2-D wave equation of the SNAPE paper, as Data.

The paper's setup:

u_tt = c (u_xx + u_yy)   on the square [-1, 1]^2,  c = 1,  t in [0, 10]
u = 0 at x = +-1 (fixed edges),  u_y = 0 at y = +-1 (free edges)
u(x, y, 0)   = 3 sin(pi x) exp(sin(pi y / 2))
u_t(x, y, 0) = atan(cos(pi x / 2))
sampled on a 50 x 50 x 100 grid

Here it is solved by expanding in the eigenfunctions of those edges, sin(m pi (x + 1) / 2) cos(n pi (y + 1) / 2), each of which oscillates at omega = (pi / 2) sqrt(c (m^2 + n^2)) -- an exact solution up to the truncation of the expansion. The paper's initial conditions are analytic, 4-periodic functions in these coordinates (odd in x + 1, even in y + 1), so their coefficients decay exponentially and 64 modes per direction reach machine precision; the generator checks the decay and refuses if the expansion has not converged.

With the default arguments the true parameters of pysnape.model.zoo.wave2d are (c, c) = (1, 1).

paper_displacement(x, y)

The paper's initial displacement, 3 sin(pi x) exp(sin(pi y / 2)).

paper_velocity(x, y)

The paper's initial velocity, atan(cos(pi x / 2)).

wave(nx=50, ny=50, nt=100, t_end=10.0, *, c=1.0, displacement=paper_displacement, velocity=paper_velocity, modes=64, quad=256)

The 2-D wave u_tt = c (u_xx + u_yy) on [-1, 1]^2, fixed at x = +-1, free at y = +-1.

Returns field u of shape (nx, ny, nt) on coordinates x, y, t. The defaults are the paper's: c = 1, a 50 x 50 x 100 grid over t in [0, 10], and its initial conditions. displacement and velocity are f(x, y) on arrays; they should satisfy the edge conditions, or the expansion converges slowly and the check refuses it.

Solver and basis

pysnape.solvers.Nested dataclass

Nested estimation settings.

Parameters:

Name	Type	Description	Default
weights	Sequence[float]	Equation weights, ascending, in balanced units (1 weighs data and equation equally). Each is one stage, warm-started from the previous one; the answer is the stage whose neighbours agree best. Too small and the data alone decide the fit; too large and the finite basis pushes the fit towards exact spline solutions of the equation, which can be trivial ones (a constant solves Burgers).	DEFAULT_WEIGHTS
starts	int	Grid starts per stage (0 disables multi-start): the starts best points of the grid, ranked by the data misfit, are searched as well as the warm start.	2
grid	Sequence[float]	Multipliers c of the starting grid theta_k = +/- c ||g_0|| / ||g_k|| at the least-squares fit, the size theta_k needs for term k to balance the parameter-free term.	(0.01, 0.1, 1.0)
max_outer	int	Levenberg-Marquardt iterations per search and the relative theta change that stops it.	40
outer_tol	int	Levenberg-Marquardt iterations per search and the relative theta change that stops it.	40
inner_tol	float	The inner fit stops when its gradient is inner_tol relative to its two parts.	1e-10
max_inner	float	The inner fit stops when its gradient is inner_tol relative to its two parts.	1e-10
stop_early	bool	Stop raising the weight once a settled answer exists (a stage whose neighbours' theta differ by less than SETTLED) and the stages after it stop being informative: two stages past it with the data misfit grown to MISFIT_GROWTH times the settled stage's and by more than MISFIT_RISE (the basis can no longer follow the data under the equation), or a stage that leaves theta exactly where it was (a stalled search). A stalled stage is always dropped from the pick. Before an answer has settled nothing stops the run: early stages may move theta far, and the misfit with it.	True
solver	str	A pysnape.solvers.linalg.METHODS entry for the Gauss-Newton matrix; "auto" is CHOLMOD when installed.	'auto'

fit_spline(asm, theta, w, beta, gtol=None, max_iter=None)

The inner problem at fixed theta: (beta, iterations).

sensitivities(asm, beta, theta, w)

d beta / d theta at the inner optimum, (n_beta, P).

Differentiating H'(H beta - y) + w J'g = 0 in theta_k gives W d beta / d theta_k = -w (J_k' g + J' g_k) with W the exact Hessian.

search(asm, theta, w, beta)

Levenberg-Marquardt on Phi_w from theta: (theta, beta, info).

run(assembly, theta0=None, beta0=None, callback=None)

Estimate the parameters from Model.unknowns (or theta0) and the least-squares fit (or beta0). callback(stage, weight, theta) is called after every stage.

pysnape.solvers.NestedResult dataclass

What one nested run produced; stage indexes the picked entry of stages.

weight property

The picked weight, in balanced units.

history property

(n_stages, n_params): theta at every weight.

pysnape.solvers.Assembly

The SNAPE system for one (model, basis, data) triple.

Attributes:

Name	Type	Description
H, y, HtH, Hty		Observation operator, observations and their normal equations, constant across sweeps. H is (n_fields * n_observed, n_beta).
n_rows		n_equations * n_points: the length of the equation residual and of the dual.
n_beta		n_fields * basis.n_basis: coefficients of all fields stacked in model order.
dtype		complex when the data or any term coefficient is complex, else float.

split(beta)

Stacked coefficient vector -> {field: its coefficients} (views).

reconstruct(beta, field, alpha=None)

Derivative alpha of field on the grid, as a grid-shaped array.

initial_beta(method='auto')

Least-squares start.

With the plain observation (H = B, no mask) this is the ordinary least-squares fit (B'B)^-1 B'y. Otherwise H'H may be singular -- the constant and linear parts of a displacement are invisible to B_tt, a masked B may leave a basis function unobserved -- so the minimum-norm least-squares solution is returned instead; the equation pins what the data cannot.

terms(beta)

(G, c) with G(beta, theta) = sum_k theta_k (G[k] @ beta + c[k]).

Index 0 is the parameter-free part; index 1 + i belongs to param_names[i]. c[k] is None when parameter k has no constant part (the common case), so the solver can skip it.

residual(beta, theta)

G(beta, theta) over all equations and grid points.

pysnape.TensorBasis

Tensor product of per-coordinate BSplineBasis objects on a fixed grid.

Parameters:

Name	Type	Description	Default
grid	Mapping[str, ArrayLike]	Ordered mapping {"x": x, "t": t} of coordinate name to its 1-D sample points. The key order is the coordinate order: the first coordinate varies fastest.	required
axes	Mapping[str, BSplineBasis]	Mapping of coordinate name to its 1-D basis; one per coordinate in grid.	required

from_spec(grid, spec) classmethod

Build the per-coordinate bases from the papers' knots / order numbers.

spec = {"x": {"knots": 200, "order": 4}, "t": {"knots": 50, "order": 3}} is the Burgers setting: knots is the number of uniformly spaced breakpoints between the smallest and largest sample of that coordinate, order the polynomial order.

alpha(alpha)

Normalise a derivative multi-index to a tuple of orders in coordinate order.

design_matrix(alpha=None)

Sparse (n_points, n_basis) matrix of derivative alpha on the grid, cached.

evaluate(points, alpha=None)

Sparse (N, n_basis) matrix of derivative alpha at scattered points.

points is {"x": xs, "t": ts} (each of length N) or an (N, ndim) array in coordinate order. Row n is the row-wise Kronecker product of the 1-D rows at point n, so on the grid points this reproduces the rows of design_matrix. Not cached.

flatten(field)

Grid array of shape grid_shape -> vector with the first coordinate fastest.

unflatten(vec)

Inverse of flatten.

coefficients(beta)

Coefficient vector -> array of shape basis_shape, same ordering as the grid.