General ODE system parameter estimator using scipy
Project description
ODE system parameter estimator
This is a little python package to provide convenience functions that use the scipy
capabilities for ODE systems and curve fitting in order to estimate the correct parameters of the ODE system for a given set of data.
Usage
In the fashion of scipy.integrate.solve_ivp
, we have an ODE system
$$
\frac{dy}{dt} = f(t,y,p)
$$
where $t$ is the independent parameter, $y$ a vector function and $p$ a parameter list. If we numerically solve the system with solve_ivp
we would get a time array solution.t
and a matrix solution.y
whose rows are the trayectories of the components of $y$.
The EstimatorProblem
class requires the data to be in this form: a time array, a samples array in row order, and the derivative function. For example, here is the Lorenz chaotic oscillator:
$$ \frac{\mathrm{d}x}{\mathrm{d}t} = \sigma (y - x), \qquad \frac{\mathrm{d}y}{\mathrm{d}t} = x (\rho - z) - y, \qquad \frac{\mathrm{d}z}{\mathrm{d}t} = x y - \beta z. $$
def lorenz(t,Y,*args):
σ,ρ,β = args
x,y,z = Y
return np.array([
σ*(y-x), # dx/dt
x*(ρ-z)-y, # dy/dt
x*y-β*z # dz/dt
])
along with some data in a file, by columns t, x, y, z
# t x y z
0.00000 1.00100 1.00100 1.00100
0.06742 1.34734 2.32878 0.95351
0.13483 2.34594 4.39342 1.17407
0.20225 4.24041 7.94361 2.22541
0.26966 7.50201 13.50709 5.74678
...
which can be imported and visualized like
t,*Y = np.loadtxt("samples.dat",unpack=True)
Y = np.array(Y)
plt.plot(t,Y.T,"o",fillstyle='none')
The parameter fitting is done by simply
problem = EstimationProblem(t,Y,lorenz)
problem.fit([9,28,2]) # initial guess
which will hopefully display a message saying that the fit was successful. That success is heavily dependant of the quality of the initial guess, specially in chaotic or ill-behaved systems.
We can test the solution with the help of the EstimationProblem.output
method, which evals the solution on the requested times:
ts = np.linspace(0,2,200)
Ys = problem.output(ts)
plt.plot(t,Y.T,"o",fillstyle='none')
plt.gca().set_prop_cycle(None) # reset color cycle
plt.plot(ts,Ys.T)
More examples
On [test/systems.py] are implementations of:
- Radioactive decay
- Housekeeping gene expression
- Lotka-Volterra
- SIR epidemics
- Damped coupled oscillators
- Lorenz-Fetter-Hamilton oscillator
All of them are tested with clean solutions and with noisy solutions to emulate experimental data.
Similar tools
There is, in Julia, a whole ecosystem of packages that implements several methods to do "dynamic data analysis". Not only there is a package named DiffEqParamEstim.jl, but there are also packages that find parts of the ODE system using scientific machine-learning.
TO DO
-
Implement dimension checking to avoid obscure error traces.
-
Finish documentation. Again.
-
Give the option to also fit the initial conditions. This should be possible, as the first data point should be a good initial guess.
-
Understand better the fail conditions, and provide more information.
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for ode_parameter_estimator-1.0.tar.gz
Algorithm | Hash digest | |
---|---|---|
SHA256 | 091bb8745cb9fb44858d9254afe4150efc4e8cc03a1b5700324c95152548803b |
|
MD5 | 28e04d8871d53402c97772960b17c720 |
|
BLAKE2b-256 | 71ef685422763b9259f36b822c898cf3724cece9074509c38e8b17a011a52929 |
Hashes for ode_parameter_estimator-1.0-py3-none-any.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 5341ad54d73a274945a05405f7c83ea5ea60eb3e688a768fb29ec3c5d08b8408 |
|
MD5 | ec90aaf8256e2579a627dfe746983b44 |
|
BLAKE2b-256 | eff9ea1873851d69edd7fc2bc8d0f3080816dc7db74ddee52d8a418959eee0f7 |