pysnape

pysnape fits the unknown parameters of an ODE or PDE directly from noisy measurements, without differentiating the noisy data and without needing initial or boundary conditions. It represents each measured field as a B-spline and estimates the spline and the model's parameters together.

The solver, Nested, works in two levels: for each candidate set of parameters it fits the spline with the equation as a penalty, and it keeps the parameters whose spline matches the data best. In our tests it stays accurate from noise-free data to heavy noise (see the examples). Statisticians know this structure as profiled estimation (Ramsay et al., 2007).

Two peer-reviewed papers back the method: MSSP 2023 (systems of ODEs and PDEs, such as Burgers, KdV, Kuramoto–Sivashinsky and coupled oscillators) and Engineering Structures 2023 (beam vibration identification with unknown boundary conditions, including a real experiment reconstructed from smartphone video).

Install

pip install pysnape

For large problems, pip install "pysnape[sparse]" adds CHOLMOD, a sparse solver that is many times faster. It needs the SuiteSparse library first (brew install suite-sparse on macOS, sudo apt install libsuitesparse-dev on Debian or Ubuntu); pysnape works without it. For plots: pip install "pysnape[plot]".

Quickstart

Write the equation the way you would on paper; pysnape reads the fields and coordinates from the data, works out which terms are unknown, chooses a basis and fits. This runs in a few seconds:

import pysnape as ps
from pysnape.data.generators import ode

# Van der Pol with mu = 2, measured with 20 % noise
data = ode.vanderpol(t_end=20.0, n=1001, mu=2.0).add_noise(percent=20, seed=2)
result = ps.fit(
    "x_tt + th1*x_t + th2*x**2*x_t + th3*x = 0",
    data,
    unknowns={"th1": 5.0, "th2": 5.0, "th3": 5.0},  # starting values
    multistart=False,  # the default, True, also tries other starts: slower, safer
)
print(result.summary())  # th1 = -2.06, th2 = 2.06, th3 = 1.02; the truth is (-2, 2, 1)

result.field() is the fitted signal on the data grid (result.field(deriv="t") its derivative), and result.plot() draws it against the measurements (needs pip install "pysnape[plot]").

How sure are the estimates? For measured data, refit with resampled residuals; for a synthetic study, refit fresh noise draws on clean data:

boot = ps.bootstrap(result, n=100, jobs=4)
print(boot.summary())  # mean, coefficient of variation and 95 % interval

clean = ps.fit(
    "x_tt + th1*x_t + th2*x**2*x_t + th3*x = 0",
    ode.vanderpol(t_end=20.0, n=1001, mu=2.0),
    unknowns={"th1": 5.0, "th2": 5.0, "th3": 5.0},
    multistart=False,
)
print(ps.ensemble(clean, n=20, noise=5, jobs=4).summary(truth={"th1": -2, "th2": 2, "th3": 1}))

The paper cases are ready-made in ps.models (Burgers, KdV, Kuramoto–Sivashinsky, Duffing, Van der Pol, the 2-D wave), as is a chain of masses, springs and dampers, ps.models.mdof_chain(n). Pass one to ps.fit in place of the equation.

The basis is chosen automatically: about one knot per three grid points, and an order of the highest derivative plus three. With the equation enforced, the basis must be fine enough to follow the solution's fastest feature. When unsure, pass candidates, knots=[100, 200, 400], and fit keeps the one that matches the data best.

Published data

The PDE-FIND datasets of Rudy et al. (Burgers, KdV, Kuramoto–Sivashinsky) load in one line:

data = ps.datasets.load("burgers")  # 256 x 101 grid
result = ps.fit(ps.models.burgers(), data)  # about two minutes
print(result.theta)  # theta1 = 0.9998, theta2 = -0.099998; the truth is (1, -0.1)
print(ps.datasets.info("burgers")["citation"])  # please cite the data's source

Large or fast-varying fields need the basis set by hand: on the chaotic Kuramoto–Sivashinsky data, knots={"x": 100, "t": 60} on a 524 x 151 window recovers the parameters within 0.5 %, where the automatic choice does not. fit warns when the fitted fields miss the data by more than 25 %, the usual sign of a basis too coarse.

Documentation

sutanubhowmick.github.io/pysnape: getting started, writing equations, data, ready-made models, how the solver works, uncertainty and the API reference. Five example notebooks run from an oscillator to published data. Every example in the docs and the notebooks is run by the test suite.

Status

Alpha: the interface may still change between minor versions. Not yet: complex-valued fields, a command-line interface with case files, and coefficients that vary in space or time. See CHANGELOG.md. Issue reports and contributions are welcome; see CONTRIBUTING.md.

Citing

If you use this software, please cite the two papers above; CITATION.cff has the details (GitHub's "Cite this repository" button reads it).

License

Apache License 2.0 — see LICENSE and NOTICE.