Profiling (ODE parameters)¶

fdatools.profiling: generalized profiling estimation of ordinary differential equation parameters (Ramsay, Hooker, Campbell & Cao 2007). ODEModel wraps a right-hand side, ProfiledODE / profile_ode fit states and parameters, and cstr_model / fitzhugh_nagumo_model are built-in examples.

Profiling: parameter estimation for nonlinear ODEs by generalized profiling.

A system of d ordinary differential equations

.. math::

\dot x_i(t) = f_i\bigl(x(t), t, \theta\bigr), \qquad i = 1, \dots, d,

is fitted to noisy observations of some (not necessarily all) of its states by the generalized profiling method of Ramsay, Hooker, Campbell & Cao (2007, JRSS-B 69: 741-796; Ramsay, Hooker & Graves 2009, chapter 11). Each state is a basis expansion :math:x_i(t) = \phi_i(t)^T c_i, and the fit has two levels.

Inner problem. For a fixed parameter vector :math:\theta the coefficients minimise the penalised criterion

.. math::

J(c \mid \theta) = \sum_{i} w_i \Bigl[ \sum_j \bigl(y_{ij} - x_i(t_{ij})\bigr)^2
    + \lambda_i \int \bigl(\dot x_i(t) - f_i(x(t), t, \theta)\bigr)^2 dt \Bigr],

where :math:w_i scales state i (R's Cwt/Twt are :math:1/w_i) and the integral is computed by a quadrature rule (composite Simpson on the B-spline break points by default, as R's quadset). J is a sum of squared residuals and is minimised by damped Gauss-Newton. An unobserved state has no data term; it is identified through the equations only.

Outer problem. The fitted coefficients :math:\hat c(\theta) are an implicit function of :math:\theta, and :math:\theta minimises the profiled data criterion

.. math::

H(\theta) = \sum_i w_i \sum_j \bigl(y_{ij} - \phi_i(t_{ij})^T \hat c_i(\theta)\bigr)^2 .

Its Jacobian comes from the implicit function theorem: at the inner optimum :math:\partial J / \partial c = 0, so

.. math::

\frac{d\hat c}{d\theta} = -\Bigl(\frac{\partial^2 J}{\partial c\,\partial c^T}\Bigr)^{-1}
    \frac{\partial^2 J}{\partial c\,\partial \theta^T},

with both second-derivative matrices exact (they include the second derivatives of f). Gauss-Newton on :math:H then gives :math:\hat\theta and its approximate covariance :math:\hat\sigma^2 (A^T A)^{-1}, :math:A = d e / d\theta, :math:\hat\sigma^2 = H(\hat\theta) / (N - p).

This replaces R's CSTR family (CSTR2in, CSTR2, CSTRfitLS, CSTRfn, CSTRres, CSTRsse) with one general interface: the user supplies the right-hand side as an :class:ODEModel. The continuously stirred tank reactor is provided as :func:cstr_model (with its input scenarios in :func:cstr_inputs), and the FitzHugh-Nagumo neuron as :func:fitzhugh_nagumo_model.

All arithmetic is done in float64 NumPy (through the array-API namespace); the estimation is an iterative optimisation, not a differentiable map, so PyTorch tensors are accepted as data but results are NumPy arrays. A right-hand side written in PyTorch is supported through :meth:ODEModel.from_torch, which obtains every derivative by automatic differentiation.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ODEModel, profile_ode
>>> decay = ODEModel(
...     rhs=lambda x, t, theta: -theta[0] * x,
...     n_states=1,
...     n_params=1,
...     jac_x=lambda x, t, theta: np.full((len(t), 1, 1), -theta[0]),
...     jac_theta=lambda x, t, theta: -x[:, :, None],
... )
>>> t = np.linspace(0.0, 2.0, 41)
>>> y = 3.0 * np.exp(-1.5 * t)
>>> basis = fdt.BSpline(domain=(0.0, 2.0), breaks=np.linspace(0.0, 2.0, 21))
>>> fit = profile_ode(decay, t, y[:, None], basis, lam=1e4, theta0=[1.0])
>>> round(float(fit.theta[0]), 4)
1.5

ODEModel `dataclass` ¶

ODEModel(
    rhs: Function,
    n_states: int,
    n_params: int,
    jac_x: Function | None = None,
    jac_theta: Function | None = None,
    hess_xx: Function | None = None,
    hess_xtheta: Function | None = None,
    state_names: tuple[str, ...] = (),
    param_names: tuple[str, ...] = (),
)

The right-hand side of an ODE system :math:\dot x = f(x, t, \theta).

Every callable receives the states x as an (n, d) array, the times t as an (n,) array and the parameters theta as a (p,) array, and is evaluated at n points at once. The derivatives are optional: any that is omitted is computed by central finite differences (of rhs for the Jacobians, of the Jacobians for the second derivatives), which is accurate to about 1e-10 relative for the Jacobians and 1e-6 for the second derivatives. Supplying them analytically, or building the model with :meth:from_torch, makes the fit exact.

Parameters:

Name	Type	Description	Default
`rhs`	`callable`	`rhs(x, t, theta)` returning `(n, d)` derivatives.	required
`n_states`	`int`	Number of states `d`.	required
`n_params`	`int`	Number of parameters `p`.	required
`jac_x`	`callable`	`(n, d, d)` array, `[q, i, k]` = :math:`\partial f_i / \partial x_k`.	`None`
`jac_theta`	`callable`	`(n, d, p)` array, `[q, i, j]` = :math:`\partial f_i / \partial \theta_j`.	`None`
`hess_xx`	`callable`	`(n, d, d, d)` array, `[q, i, k, l]` = :math:`\partial^2 f_i / \partial x_k \partial x_l`.	`None`
`hess_xtheta`	`callable`	`(n, d, d, p)` array, `[q, i, k, j]` = :math:`\partial^2 f_i / \partial x_k \partial \theta_j`.	`None`
`state_names`	`sequence of str`	Labels; default `x0, x1, ...` and `theta0, theta1, ...`.	`()`
`param_names`	`sequence of str`	Labels; default `x0, x1, ...` and `theta0, theta1, ...`.	`()`

Raises:

Type	Description
`ValueError`	If a count is not positive or a name list has the wrong length.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import ODEModel
>>> growth = ODEModel(lambda x, t, th: th[0] * x, n_states=1, n_params=1)
>>> growth(np.array([[2.0]]), np.array([0.0]), np.array([0.5])).tolist()
[[1.0]]
>>> round(
...     float(
...         growth.jacobians(np.array([[2.0]]), np.array([0.0]), np.array([0.5]))[1][0, 0, 0]
...     ),
...     6,
... )
2.0

call ¶

__call__(x: Any, t: Any, theta: Any) -> Array

Evaluate :math:f(x, t, \theta) at n points.

Parameters:

Name	Type	Description	Default
`x`	`array_like`	States, `(n, d)`.	required
`t`	`array_like`	Times, `(n,)`.	required
`theta`	`array_like`	Parameters, `(p,)`.	required

Returns:

Type	Description
`array`	`(n, d)` derivatives.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import fitzhugh_nagumo_model
>>> fhn = fitzhugh_nagumo_model()
>>> fhn(np.array([[0.0, 0.0]]), np.array([0.0]), np.array([0.2, 0.2, 3.0])).round(
...     4
... ).tolist()
[[0.0, 0.0667]]

jacobians ¶

jacobians(
    x: Any, t: Any, theta: Any
) -> tuple[Array, Array]

Return :math:\partial f / \partial x and :math:\partial f / \partial \theta.

Parameters:

Name	Type	Description	Default
`x`	`array_like`	States, `(n, d)`.	required
`t`	`array_like`	Times, `(n,)`.	required
`theta`	`array_like`	Parameters, `(p,)`.	required

Returns:

Name	Type	Description
`jac_x`	`array`	`(n, d, d)`.
`jac_theta`	`array`	`(n, d, p)`.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import fitzhugh_nagumo_model
>>> fx, ft = fitzhugh_nagumo_model().jacobians(
...     np.array([[1.0, 0.5]]), np.array([0.0]), np.array([0.2, 0.2, 3.0])
... )
>>> fx.shape, ft.shape
((1, 2, 2), (1, 2, 3))

hessians ¶

hessians(x: Any, t: Any, theta: Any) -> tuple[Array, Array]

Return the second derivatives of f used by the profiling fit.

Parameters:

Name	Type	Description	Default
`x`	`array_like`	States, `(n, d)`.	required
`t`	`array_like`	Times, `(n,)`.	required
`theta`	`array_like`	Parameters, `(p,)`.	required

Returns:

Name	Type	Description
`hess_xx`	`array`	`(n, d, d, d)`, :math:`\partial^2 f_i / \partial x_k \partial x_l`.
`hess_xtheta`	`array`	`(n, d, d, p)`, :math:`\partial^2 f_i / \partial x_k \partial \theta_j`.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import fitzhugh_nagumo_model
>>> fxx, fxt = fitzhugh_nagumo_model().hessians(
...     np.array([[1.0, 0.5]]), np.array([0.0]), np.array([0.2, 0.2, 3.0])
... )
>>> float(fxx[0, 0, 0, 0])
-6.0

simulate ¶

simulate(
    t: Any,
    x0: Any,
    theta: Any,
    *,
    breaks: Sequence[float] = (),
    rtol: float = _ODE_RTOL,
    atol: float = _ODE_ATOL,
    method: OdeMethod = "LSODA",
) -> Array

Solve the initial-value problem with :func:scipy.integrate.solve_ivp.

Parameters:

Name	Type	Description	Default
`t`	`array_like`	Increasing output times; the integration starts at `t[0]`.	required
`x0`	`array_like`	`(d,)` state at `t[0]`.	required
`theta`	`array_like`	`(p,)` parameters.	required
`breaks`	`sequence of float`	Times where the right-hand side jumps (for example step inputs). The integration is restarted at each, so the solver never steps across a discontinuity.	`()`
`rtol`	`float`	Solver tolerances.	`_ODE_RTOL`
`atol`	`float`	Solver tolerances.	`_ODE_RTOL`
`method`	`str`	A :func:`~scipy.integrate.solve_ivp` method; LSODA by default.	`'LSODA'`

Returns:

Type	Description
`array`	`(len(t), d)` solution.

Raises:

Type	Description
`ValueError`	If `t` is not increasing, a shape is wrong, or the solver fails.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import ODEModel
>>> decay = ODEModel(lambda x, t, th: -th[0] * x, n_states=1, n_params=1)
>>> x = decay.simulate([0.0, 1.0], [1.0], [2.0])
>>> round(float(x[1, 0]), 8) == round(float(np.exp(-2.0)), 8)
True

from_torch `classmethod` ¶

from_torch(
    rhs: Callable[[Any, Any, Any], Any],
    n_states: int,
    n_params: int,
    *,
    state_names: Sequence[str] = (),
    param_names: Sequence[str] = (),
) -> ODEModel

Build a model from a PyTorch right-hand side, with autodiff derivatives.

rhs(x, t, theta) must use PyTorch operations only and treat the n points independently (it is differentiated one point at a time with :func:torch.func.vmap). Every Jacobian and second derivative is then exact. Requires the fdatools[torch] extra; PyTorch is imported only when this method is called.

Parameters:

Name	Type	Description	Default
`rhs`	`callable`	PyTorch function of `(x, t, theta)` returning `(n, d)` tensors.	required
`n_states`	`int`	Dimensions of `x` and `theta`.	required
`n_params`	`int`	Dimensions of `x` and `theta`.	required
`state_names`	`sequence of str`	Labels.	`()`
`param_names`	`sequence of str`	Labels.	`()`

Returns:

Type	Description
`ODEModel`	A model whose callables accept and return NumPy arrays.

Examples:

>>> import numpy as np
>>> import torch
>>> from fdatools.profiling import ODEModel
>>> model = ODEModel.from_torch(lambda x, t, th: -th[0] * x**2, 1, 1)
>>> float(
...     model.hessians(np.array([[3.0]]), np.array([0.0]), np.array([2.0]))[0][0, 0, 0, 0]
... )
-4.0

Residuals ¶

Bases: NamedTuple

Residuals of the inner criterion and their Jacobians at (c, theta).

J(c | theta) is data @ data + equation @ equation. Rows of the Jacobians are the data residuals (observed states in order) followed by the equation residuals (state by state, quadrature node by node); columns of jac_coefs are the stacked coefficients of all states.

Attributes:

Name	Type	Description
`data`	`array`	`sqrt(w_i) (y_ij - x_i(t_ij))`.
`equation`	`array`	`sqrt(w_i lambda_i v_q) (Dx_i(s_q) - f_i(x(s_q), s_q, theta))`.
`jac_coefs`	`array`	Derivative of `concat(data, equation)` with respect to the coefficients.
`jac_theta`	`array`	Derivative of `concat(data, equation)` with respect to `theta`.

InnerFit `dataclass` ¶

InnerFit(
    coefs: tuple[Array, ...],
    criterion: float,
    sse: float,
    df: float,
    gcv: float,
    dcoefs_dtheta: Array,
    n_iter: int,
    converged: bool,
)

The state fit for one fixed parameter vector (the inner problem).

Attributes:

Name	Type	Description
`coefs`	`tuple of array`	Fitted basis coefficients, one `(n_basis_i,)` array per state.
`criterion`	`float`	The minimised `J(c \| theta)`.
`sse`	`float`	Its weighted data part `sum_i w_i sum_j (y_ij - x_i(t_ij))**2`.
`df`	`float`	Effective degrees of freedom, the trace of the (linearised) map from the data to the fitted values.
`gcv`	`float`	`N * sse / (N - df)**2` for `N` observations.
`dcoefs_dtheta`	`array`	`(K, p)` implicit-function derivative of the stacked coefficients.
`n_iter`	`int`	Gauss-Newton iterations used.
`converged`	`bool`	Whether the step tolerance was met.

ProfileResult `dataclass` ¶

ProfileResult(
    theta: Array,
    cov: Array,
    sigma2: float,
    states: tuple[FData, ...],
    inner: InnerFit,
    theta_path: Array,
    sse_path: Array,
    n_iter: int,
    converged: bool,
    param_names: tuple[str, ...] = (),
    state_names: tuple[str, ...] = (),
)

The result of a generalized-profiling fit.

Attributes:

Name	Type	Description
`theta`	`array`	`(p,)` estimated parameters.
`cov`	`array`	`(p, p)` approximate covariance `sigma2 * inv(A.T @ A)`.
`sigma2`	`float`	Residual variance `sse / (N - p)`.
`states`	`tuple of FData`	The fitted state functions, one single-curve :class:`~fdatools.FData` per state.
`inner`	`InnerFit`	The inner fit at `theta`.
`theta_path`	`array`	`(n_iter + 1, p)` parameter values visited, from the start.
`sse_path`	`array`	`(n_iter + 1,)` profiled data criterion at each of them.
`n_iter`	`int`	Outer Gauss-Newton iterations.
`converged`	`bool`	Whether the outer tolerance was met.
`param_names, state_names`	`tuple of str`	Labels from the model.

stderr `property` ¶

stderr: Array

Standard errors, the square roots of the diagonal of :attr:cov.

Returns:

Type	Description
`array`	`(p,)` standard errors.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ODEModel, profile_ode
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 21)
>>> y = np.exp(-t) + 0.01 * np.cos(40 * t)
>>> fit = profile_ode(m, t, y[:, None], fdt.BSpline((0.0, 1.0), 12), 1e3, [0.5])
>>> fit.stderr.shape
(1,)

call ¶

__call__(t: Any, deriv: int = 0) -> Array

Evaluate the fitted states at t.

Parameters:

Name	Type	Description	Default
`t`	`array_like`	Evaluation times.	required
`deriv`	`int`	Derivative order.	`0`

Returns:

Type	Description
`array`	`(len(t), d)` state values.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ODEModel, profile_ode
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 21)
>>> fit = profile_ode(m, t, np.exp(-t)[:, None], fdt.BSpline((0.0, 1.0), 12), 1e3, [0.5])
>>> fit(np.array([0.0, 1.0])).shape
(2, 1)

ProfiledODE ¶

ProfiledODE(
    model: ODEModel,
    t: Any,
    y: Any,
    bases: Basis | Sequence[Basis],
    lam: float | Sequence[float] = 1.0,
    *,
    state_weights: float | Sequence[float] = 1.0,
    quadrature: tuple[Any, Any] | None = None,
    n_quad: int = 5,
)

A generalized-profiling problem: model, data, bases and smoothing.

Parameters:

Name	Type	Description	Default
`model`	`ODEModel`	The right-hand side.	required
`t`	`array_like or sequence`	Observation times: one `(n,)` array shared by all states, or one array per state (`None` for an unobserved state).	required
`y`	`array_like or sequence`	Observations: an `(n, d)` array (`NaN` marks a missing value; an all-`NaN` column an unobserved state), or one `(n_i,)` array per state (`None` for an unobserved state).	required
`bases`	`Basis or sequence of Basis`	One basis for every state, or one per state, all on the same domain.	required
`lam`	`float or sequence of float`	Weight :math:`\lambda_i` of the equation term of each state.	`1.0`
`state_weights`	`float or sequence of float`	Weights :math:`w_i` multiplying both terms of each state (the reciprocal of R's `Cwt`/`Twt`); `1` by default. Use them to put states measured on different scales on a common footing.	`1.0`
`quadrature`	`tuple of array_like`	`(nodes, weights)` for the integral. By default the composite Simpson rule (:func:`simpson_rule`) with `n_quad` nodes on each panel between the B-spline break points of the bases (R's `quadset`).	`None`
`n_quad`	`int`	Nodes per panel of the default rule.	`5`

Attributes:

Name	Type	Description
`model`	`ODEModel`
`bases`	`tuple of Basis`
`lam, state_weights`	`tuple of float`
`nodes, weights`	`array`	The quadrature rule.
`n_obs`	`int`	Total number of observations `N`.

Raises:

Type	Description
`ValueError`	If shapes, domains or weights are inconsistent, or no state is observed.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, fitzhugh_nagumo_model
>>> fhn = fitzhugh_nagumo_model()
>>> t = np.linspace(0.0, 20.0, 201)
>>> x = fhn.simulate(t, [-1.0, 1.0], [0.2, 0.2, 3.0])
>>> basis = fdt.BSpline(domain=(0.0, 20.0), breaks=np.linspace(0.0, 20.0, 201))
>>> problem = ProfiledODE(fhn, t, x, basis, lam=1e3)
>>> fit = problem.fit([0.3, 0.3, 2.5])
>>> np.round(fit.theta, 3).tolist()
[0.2, 0.2, 3.0]

n_coefs `property` ¶

n_coefs: int

Total number of basis coefficients K over all states.

Returns:

Type	Description
`int`	`sum(b.n_basis for b in bases)`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, fitzhugh_nagumo_model
>>> t = np.linspace(0.0, 1.0, 5)
>>> p = ProfiledODE(fitzhugh_nagumo_model(), t, np.ones((5, 2)), fdt.BSpline(n_basis=6))
>>> p.n_coefs
12

residuals ¶

residuals(coefs: Any, theta: Any) -> Residuals

Residuals and Jacobians of the inner criterion (R's CSTRfitLS).

Parameters:

Name	Type	Description	Default
`coefs`	`array_like or sequence of array_like`	The stacked coefficients, or one array per state.	required
`theta`	`array_like`	`(p,)` parameters.	required

Returns:

Type	Description
`Residuals`	`data`, `equation`, `jac_coefs`, `jac_theta`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, ODEModel
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 5)
>>> p = ProfiledODE(m, t, np.ones((5, 1)), fdt.BSpline(n_basis=4))
>>> r = p.residuals(np.ones(4), [0.0])
>>> float(abs(r.data).max()), float(abs(r.equation).max())
(0.0, 0.0)

criterion ¶

criterion(coefs: Any, theta: Any) -> float

Evaluate the inner criterion J(c | theta).

Parameters:

Name	Type	Description	Default
`coefs`	`array_like or sequence of array_like`	Coefficients.	required
`theta`	`array_like`	Parameters.	required

Returns:

Type	Description
`float`	Sum of squared data and equation residuals.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, ODEModel
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 5)
>>> p = ProfiledODE(m, t, np.zeros((5, 1)), fdt.BSpline(n_basis=4))
>>> p.criterion(np.ones(4), [0.0])
5.0

fit_states ¶

fit_states(
    theta: Any,
    coef0: Any = None,
    *,
    tol: float = 1e-14,
    max_iter: int = 200,
) -> InnerFit

Fit the states for fixed parameters (the inner problem, R's CSTRfn).

Parameters:

Name	Type	Description	Default
`theta`	`array_like`	`(p,)` parameters.	required
`coef0`	`array_like or sequence`	Starting coefficients; by default a least-squares fit of each observed state and zeros for unobserved ones.	`None`
`tol`	`float`	Step size, relative to each state's largest coefficient, at which Gauss-Newton stops.	`1e-14`
`max_iter`	`int`	Iteration cap.	`200`

Returns:

Type	Description
`InnerFit`	Coefficients, criterion, degrees of freedom and `dc/dtheta`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, ODEModel
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 11)
>>> p = ProfiledODE(m, t, np.exp(-t)[:, None], fdt.BSpline(n_basis=8), lam=1e2)
>>> inner = p.fit_states([1.0])
>>> inner.converged, inner.sse < 1e-8
(True, True)

fit ¶

fit(
    theta0: Any,
    coef0: Any = None,
    *,
    tol: float = 1e-08,
    max_iter: int = 200,
    inner_tol: float = 1e-14,
    inner_max_iter: int = 200,
) -> ProfileResult

Estimate the parameters by Gauss-Newton on the profiled criterion.

Each outer step re-fits the states (warm-started from the previous fit), forms the data-residual Jacobian A = -Phi dc/dtheta by the implicit function theorem, and takes a Gauss-Newton step on theta with Marquardt damping when the criterion would increase. Iteration stops when the relative-offset criterion of Bates and Watts (as R's nls), |Q'e| / sqrt(p) / (|e| / sqrt(N - p)), falls below tol, or the step is at rounding level.

Parameters:

Name	Type	Description	Default
`theta0`	`array_like`	`(p,)` starting parameters.	required
`coef0`	`array_like or sequence`	Starting coefficients for the first inner fit.	`None`
`tol`	`float`	Relative-offset tolerance.	`1e-08`
`max_iter`	`int`	Outer iteration cap.	`200`
`inner_tol`	`optional`	Passed to the inner fits (see :meth:`fit_states`).	`1e-14`
`inner_max_iter`	`optional`	Passed to the inner fits (see :meth:`fit_states`).	`1e-14`

Returns:

Type	Description
`ProfileResult`	Estimates, covariance, fitted states and the iteration path.

Raises:

Type	Description
`ValueError`	If there are no more observations than parameters.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import ProfiledODE, ODEModel
>>> m = ODEModel(lambda x, t, th: -th[0] * x, 1, 1)
>>> t = np.linspace(0.0, 1.0, 21)
>>> basis = fdt.BSpline(domain=(0.0, 1.0), breaks=np.linspace(0.0, 1.0, 11))
>>> fit = ProfiledODE(m, t, 2 * np.exp(-0.7 * t)[:, None], basis, lam=1e4).fit([0.2])
>>> round(float(fit.theta[0]), 5)
0.7

fitzhugh_nagumo_model ¶

fitzhugh_nagumo_model() -> ODEModel

Build the FitzHugh-Nagumo neuron model, with analytic derivatives.

.. math::

\dot V = c\,(V - V^3/3 + R), \qquad \dot R = -(V - a + b R) / c,

with parameters :math:\theta = (a, b, c); the classic values are (0.2, 0.2, 3) (Ramsay et al. 2007, section 2.1).

Returns:

Type	Description
`ODEModel`	States `("V", "R")`, parameters `("a", "b", "c")`.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import fitzhugh_nagumo_model
>>> fhn = fitzhugh_nagumo_model()
>>> x = fhn.simulate(np.linspace(0.0, 20.0, 5), [-1.0, 1.0], [0.2, 0.2, 3.0])
>>> x.shape
(5, 2)

cstr_inputs ¶

cstr_inputs(
    t: Any, condition: str = "all.cool.step"
) -> Array

Input series of the CSTR experiments (R's CSTR2in).

The reactor is driven by five inputs: the flow F, the inlet concentration CA0, the inlet temperature T0, the coolant inlet temperature Tcin and the coolant flow Fc. The all.*.step scenarios hold them at nominal values (F = 1, CA0 = 2, T0 = 323, Fc = 15, Tcin = 335 cool or 365 hot) and step one input at a time on the 4-minute intervals [4k, 4k + 4) of [0, 64): F to 1.5 and 0.5, CA0 to 2.2 and 1.8, T0 to 343 and 303, Tcin up and down by 5, Fc to 20 and 10. The Tc.*.step scenarios lower Tcin by 15 on [2, 12) only, with F = 1 (hot) or 0.05 (cool).

Parameters:

Name	Type	Description	Default
`t`	`array_like`	Times (minutes).	required
`condition`	`str`	One of :data:`CSTR_CONDITIONS`.	`'all.cool.step'`

Returns:

Type	Description
`array`	`(len(t), 5)` inputs, columns :data:`CSTR_INPUTS`.

Raises:

Type	Description
`ValueError`	If `condition` is unknown.

Examples:

>>> from fdatools.profiling import cstr_inputs
>>> cstr_inputs([0.0, 4.0, 8.0], "all.cool.step")[:, 0].tolist()
[1.0, 1.5, 0.5]

cstr_model ¶

cstr_model(
    condition: str
    | Callable[[Array], Any] = "all.cool.step",
    *,
    estimate: Sequence[str] = ("kref", "EoverR", "a", "b"),
    **values: float,
) -> ODEModel

Build the two-state continuously stirred tank reactor (Ramsay et al. 2007).

States are the concentration C and temperature T:

.. math::

\dot C &= -\bigl(k(T) + F/V\bigr) C + (F/V)\, C_{A0}, \\
\dot T &= -\bigl(\alpha + F/V\bigr) T - \frac{\Delta H}{V \rho C_p} k(T)\, C
          + (F/V)\, T_0 + \alpha\, T_{cin},

with :math:k(T) = k_{ref} \exp\{-10^4 (E/R) (1/T - 1/T_{ref})\} and :math:\alpha = a F_c^{b+1} / \{V \rho C_p (F_c + a F_c^b / (2 \rho_c C_{pc}))\} -- the same right-hand side as R's CSTR2. The parameters are kref, EoverR (in units of 1e4 K), a and b; those named in estimate form :math:\theta (in that order) and the others are held at their values. All derivatives are analytic.

Parameters:

Name	Type	Description	Default
`condition`	`str or callable`	A scenario of :func:`cstr_inputs`, or a function mapping `(n,)` times to the `(n, 5)` inputs (columns :data:`CSTR_INPUTS`).	`'all.cool.step'`
`estimate`	`sequence of str`	The parameters to estimate.	`('kref', 'EoverR', 'a', 'b')`
`**values`	`float`	Overrides of the parameter values (:data:`CSTR_PARAMETERS`, used for the ones not estimated) and of the constants (:data:`CSTR_CONSTANTS`).	`{}`

Returns:

Type	Description
`ODEModel`	States `("C", "T")`, parameters `estimate`.

Raises:

Type	Description
`ValueError`	If a name in `estimate` or `values` is unknown, or `estimate` is empty or repeats a name.

Examples:

>>> import numpy as np
>>> from fdatools.profiling import cstr_model
>>> model = cstr_model("all.cool.step", estimate=("kref", "EoverR"))
>>> dx = model(np.array([[1.5965, 341.3754]]), np.array([1.0]), np.array([0.461, 0.83301]))
>>> dx.round(6).tolist()
[[0.000106, 0.010411]]

simpson_rule ¶

simpson_rule(
    breaks: Sequence[float], n_quad: int = 5
) -> tuple[Array, Array]

Composite Simpson rule on the panels between consecutive breaks.

Each panel gets n_quad equally spaced nodes including both ends, with weights h/3 * (1, 4, 2, 4, ..., 2, 4, 1); a node on an interior break appears twice, once per panel -- the rule of R's quadset.

Parameters:

Name	Type	Description	Default
`breaks`	`sequence of float`	Strictly increasing panel boundaries.	required
`n_quad`	`int`	Odd number of nodes per panel, at least 3.	`5`

Returns:

Type	Description
`nodes, weights : array`	`((len(breaks) - 1) * n_quad,)` each.

Raises:

Type	Description
`ValueError`	If `n_quad` is even or below 3, or `breaks` is not increasing.

Examples:

>>> from fdatools.profiling import simpson_rule
>>> nodes, weights = simpson_rule([0.0, 1.0], 5)
>>> nodes.tolist()
[0.0, 0.25, 0.5, 0.75, 1.0]
>>> [round(float(w) * 12, 12) for w in weights]
[1.0, 4.0, 2.0, 4.0, 1.0]

profile_ode ¶

profile_ode(
    model: ODEModel,
    t: Any,
    y: Any,
    bases: Basis | Sequence[Basis],
    lam: float | Sequence[float],
    theta0: Any,
    *,
    state_weights: float | Sequence[float] = 1.0,
    quadrature: tuple[Any, Any] | None = None,
    n_quad: int = 5,
    coef0: Any = None,
    tol: float = 1e-08,
    max_iter: int = 200,
) -> ProfileResult

Estimate ODE parameters from noisy data by generalized profiling.

A one-call front end to :class:ProfiledODE; see there and the module documentation for the method and the arguments.

Parameters:

Name	Type	Description	Default
`model`	`ODEModel`	The right-hand side `f(x, t, theta)`.	required
`t`	`array_like or sequence`	Observations (see :class:`ProfiledODE`).	required
`y`	`array_like or sequence`	Observations (see :class:`ProfiledODE`).	required
`bases`	`Basis or sequence of Basis`	State bases.	required
`lam`	`float or sequence of float`	Equation weights.	required
`theta0`	`array_like`	Starting parameters.	required
`state_weights`	`float or sequence of float`	Per-state weights.	`1.0`
`quadrature`	`tuple of array_like`	`(nodes, weights)`.	`None`
`n_quad`	`int`	Simpson nodes per panel.	`5`
`coef0`	`array_like or sequence`	Starting coefficients.	`None`
`tol`	`float`	Outer relative-offset tolerance.	`1e-08`
`max_iter`	`int`	Outer iteration cap.	`200`

Returns:

Type	Description
`ProfileResult`	Estimates, covariance, fitted states and the iteration path.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.profiling import fitzhugh_nagumo_model, profile_ode
>>> fhn = fitzhugh_nagumo_model()
>>> t = np.linspace(0.0, 20.0, 201)
>>> x = fhn.simulate(t, [-1.0, 1.0], [0.2, 0.2, 3.0])
>>> y = np.column_stack([x[:, 0], np.full(201, np.nan)])  # R is never observed
>>> basis = fdt.BSpline(domain=(0.0, 20.0), breaks=np.linspace(0.0, 20.0, 201))
>>> fit = profile_ode(fhn, t, y, basis, lam=1e3, theta0=[0.25, 0.25, 2.8])
>>> np.round(fit.theta, 2).tolist()
[0.2, 0.2, 3.0]

Profiling (ODE parameters)¶

ODEModel dataclass ¶

__call__ ¶

jacobians ¶

hessians ¶

simulate ¶

from_torch classmethod ¶

Residuals ¶

InnerFit dataclass ¶

ProfileResult dataclass ¶

stderr property ¶

__call__ ¶

ProfiledODE ¶

n_coefs property ¶

residuals ¶

criterion ¶

fit_states ¶

fit ¶

fitzhugh_nagumo_model ¶

cstr_inputs ¶

cstr_model ¶

simpson_rule ¶

profile_ode ¶

ODEModel `dataclass` ¶

call ¶

from_torch `classmethod` ¶

InnerFit `dataclass` ¶

ProfileResult `dataclass` ¶

stderr `property` ¶

call ¶

n_coefs `property` ¶