Dynamics¶

fdatools.dynamics: principal differential analysis with the PDA estimator (fit, ODE solve, overlay plot) and the phase_plane() plot.

Dynamics: principal differential analysis and phase-plane plots.

Principal differential analysis (PDA, Ramsay & Silverman chapter 19; Ramsay, Hooker & Graves chapter 11) looks for the linear differential equation that a sample of curves satisfies as closely as possible. For a single variable and order m it estimates weight functions :math:\beta_0, \dots, \beta_{m-1} such that

.. math::

L x_n = D^m x_n + \sum_{j=0}^{m-1} \beta_j(t)\, D^j x_n \approx 0
\qquad (n = 1, \dots, N),

by minimising

.. math::

\frac{1}{N} \sum_{n=1}^{N} \int (L x_n)(t)^2 \, dt
+ \sum_{j} \lambda_j \int \bigl(K \beta_j\bigr)(t)^2 \, dt ,

where :math:K is a roughness operator on the weights. Expanding each weight in a basis, :math:\beta_j = \theta_j^{T} b_j, makes this a linear least-squares problem in the stacked coefficients :math:b:

.. math::

\Bigl(A + \operatorname{diag}_j \lambda_j R_j\Bigr) b = r , \quad
A_{(j,a),(l,c)} = \frac{1}{N} \sum_n \int \theta_{ja}\theta_{lc}\,
    D^j x_n\, D^l x_n , \quad
r_{(j,a)} = -\frac{1}{N} \sum_n \int \theta_{ja}\, D^j x_n\, D^m x_n .

A system of d coupled variables is the same problem once per equation i, with the unknowns :math:\beta_{ikj} multiplying :math:D^j x_k for every variable k.

Forcing functions :math:u_1, \dots, u_K (inputs from outside the system) enter an equation through weight functions :math:\alpha_1, \dots, \alpha_K,

.. math::

D^m x = -\sum_{j} \beta_j D^j x + \sum_{k} \alpha_k u_k ,

so the residual becomes :math:L x - \sum_k \alpha_k u_k and each :math:\alpha_k u_k is one more regressor, :math:-u_k, in the same least-squares problem (fitted jointly with the :math:\beta_j).

The fitted equation, written as a first-order system :math:Dz = A(t) z + f(t) in the state :math:z = (x, Dx, \dots, D^{m-1} x), has a local stability picture: the eigenvalues of :math:A(t) (negative real parts decay, imaginary parts oscillate) and the equilibrium :math:z^* = -A(t)^{-1} f(t) the system is drawn towards (:meth:PDA.stability).

Replaces R's pda.fd (:class:PDA, including awtlist/ufdlist), eigen.pda (:meth:PDA.stability), pda.overlay (:meth:PDA.plot_overlay) and phaseplanePlot (:func:phase_plane); the fitted equation is integrated with :func:scipy.integrate.solve_ivp (:meth:PDA.solve).

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> curves = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.5], [0.0, 2.0]]), basis)
>>> pda = PDA(order=2, n_grid=None).fit(curves)  # a sin t + b cos t solve D²x + x = 0
>>> [round(float(w.coefs[0, 0]), 10) + 0.0 for w in pda.weights_]
[1.0, 0.0]

PDAStability `dataclass` ¶

PDAStability(t: Array, eigenvalues: Array, limits: Array)

Local stability of a fitted differential equation over time.

Replaces the value of R's eigen.pda. At each time t the fitted equations are written as the first-order system Dz = A(t) z + f(t) in the state z = (x_1, Dx_1, …, D^{m-1}x_1, x_2, …), where A(t) is the companion matrix of the weights β and f(t) holds the forcing terms Σ_k a_k(t) u_k(t).

Attributes:

Name	Type	Description
`t`	`array`	The `(n_t,)` time points.
`eigenvalues`	`array`	`(n_t, n_vars * order)` complex eigenvalues of `A(t)`, each row in decreasing modulus (a complex pair keeps its positive imaginary part first). A negative real part means local exponential decay, a positive one growth, and a non-zero imaginary part oscillation with angular frequency `\|Im λ\|`.
`limits`	`array`	`(n_t, n_vars * order)` equilibrium states `z* = -A(t)⁻¹ f(t)`, the point the system is drawn to (when it is stable) if the weights and forcing were frozen at their value at `t`. Zero for an unforced equation and `nan` where `A(t)` is exactly singular.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> result = PDA(order=2, n_grid=None).fit(fd).stability(n_points=3)
>>> np.round(result.eigenvalues[0], 10) + 0.0
array([0.+1.j, 0.-1.j])

plot ¶

plot(ax: Axes | None = None, **kwargs: Any) -> Axes

Plot the real (solid) and imaginary (dashed) parts of the eigenvalues.

Parameters:

Name	Type	Description	Default
`ax`	`Axes`	Axes to draw on. A new figure is created when omitted.	`None`
`**kwargs`	`Any`	Passed to :meth:`matplotlib.axes.Axes.plot` for every line.	`{}`

Returns:

Type	Description
`Axes`	The axes drawn on: one solid line per eigenvalue for its real part, then one dashed line per eigenvalue for its imaginary part.

Examples:

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> len(PDA(order=2).fit(fd).stability().plot().lines)
5

PDA ¶

PDA(
    order: int = 2,
    *,
    weight_basis: Basis | Sequence[Basis] | None = None,
    lam: float | Sequence[float] = 0.0,
    penalty: int | LDO = 2,
    n_grid: int | None = _R_GRID,
    forcing_basis: Basis | Sequence[Any] | None = None,
    forcing_lam: float | Sequence[Any] = 0.0,
)

Bases: TransformerMixin, BaseEstimator

Principal differential analysis: fit a linear ODE to a sample of curves.

Replaces R's pda.fd. Estimates weight functions β_j so that D^m x + Σ_j β_j D^j x is as small as possible over the sample, in the mean integrated square, with an optional roughness penalty on each weight (see the module docstring for the criterion). Curves with several variables (n_vars = d > 1) define a system of d coupled equations. Forcing functions u_k passed to :meth:fit add terms a_k u_k to the right-hand side, D^m x = -Σ_j β_j D^j x + Σ_k a_k u_k (R's awtlist/ufdlist); their weights a_k are estimated jointly with the β_j.

Parameters:

Name	Type	Description	Default
`order`	`int`	Order `m` of the equation. Default `2`.	`2`
`weight_basis`	`Basis or sequence of Basis`	Basis for the weight functions; one per derivative order `j` when a sequence of length `order`. `None` (default) uses a constant basis, i.e. a constant-coefficient equation.	`None`
`lam`	`float or sequence of float`	Roughness penalty on each weight function, per derivative order when a sequence. Default `0.0`.	`0.0`
`penalty`	`int or LDO`	Roughness operator applied to the weights (`β` and `a`). Default `2`.	`2`
`n_grid`	`int or None`	How the integrals are computed. An integer is the number of equally spaced points of a trapezoidal rule; the default `501` is the rule R's `pda.fd` uses, so fdatools reproduces it. (R raises its grid to five times the number of basis functions of the curves when that is larger; pass that number to reproduce R for curve bases of more than 100 functions.) The residual functions are then least-squares fits of the residual on the same grid. `None` integrates exactly (composite Gauss-Legendre on the break points, exact for spline curves and weights) and projects the residuals in `L²`; it differs from R's rule by about `1e-5` relative on realistic data.	`_R_GRID`
`forcing_basis`	`Basis or sequence`	Basis for the forcing weights `a`. `None` (default) is a constant basis. A single basis applies to every forcing function; otherwise one entry per forcing function (for a system: one entry per equation, each a basis or one basis per forcing function of that equation).	`None`
`forcing_lam`	`float or sequence`	Roughness penalty on each forcing weight, spread like `forcing_basis`. Default `0.0`.	`0.0`

Attributes:

Name	Type	Description
`weights_`	`tuple`	For one variable, `weights_[j]` is `β_j` as a one-curve :class:`~fdatools.core.FData`. For `d` variables, `weights_[i][k][j]` multiplies `D^j x_k` in equation `i`. This mirrors the nesting of R's `bwtlist`.
`forcing_weights_`	`tuple`	For one variable, `forcing_weights_[k]` is `a_k`; for `d` variables, `forcing_weights_[i][k]` multiplies forcing function `k` of equation `i` (R's `awtlist`). Empty when fitted without forcing.
`residuals_`	`FData`	`L x - Σ_k a_k u_k` for the fitted curves, in their basis.
`operator_`	`LDO or None`	The fitted operator `L` for one variable (so that `fd(t, pda.operator_)` evaluates `L x`); `None` for a system.
`n_vars_`	`int`	Number of variables of the fitted curves.
`domain_`	`tuple of float`	Domain of the fitted curves.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> t = np.linspace(0.0, 1.0, 101)
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=24, order=5)
>>> from fdatools.smoothing import smooth
>>> fd = smooth(np.exp(-4 * t), t, basis=basis, lam=0.0).fd
>>> pda = PDA(order=1).fit(fd)  # Dx + βx = 0
>>> round(float(pda.weights_[0].coefs[0, 0]), 6)
4.0

A forcing function: x = 0.5 (1 - exp(-4t)) solves Dx = -4x + 2u for the constant input u = 1.

>>> fd = smooth(0.5 * (1 - np.exp(-4 * t)), t, basis=basis, lam=0.0).fd
>>> u = fdt.FData(np.array([1.0]), fdt.Constant(domain=(0.0, 1.0)))
>>> pda = PDA(order=1).fit(fd, forcing=u)
>>> (
...     round(float(pda.weights_[0].coefs[0, 0]), 6),
...     round(float(pda.forcing_weights_[0].coefs[0, 0]), 6),
... )
(4.0, 2.0)

fit ¶

fit(X: FData, y: Any = None, *, forcing: Any = None) -> PDA

Estimate the weight functions from a sample of curves.

Parameters:

Name	Type	Description	Default
`X`	`FData`	The curves; a multivariate :class:`~fdatools.core.FData` (`n_vars > 1`) is fitted as a system of coupled equations.	required
`y`	`None`	Ignored; present for scikit-learn compatibility.	`None`
`forcing`	`FData or sequence`	Forcing functions `u_k` (R's `ufdlist`). For one equation an :class:`~fdatools.core.FData` or a sequence of them; for a system one entry per equation, each `None`, an :class:`~fdatools.core.FData` or a sequence of them. Each has one variable and either one curve per curve of `X` or a single curve shared by all of them.	`None`

Returns:

Type	Description
`PDA`	The fitted estimator.

Raises:

Type	Description
`TypeError`	If `X` or a forcing function is not an :class:`~fdatools.core.FData`.
`ValueError`	If a parameter is invalid or the normal equations are singular.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> len(PDA(order=2).fit(fd).weights_)
2

transform ¶

transform(X: FData, *, forcing: Any = None) -> FData

Apply the fitted equation to curves: return L x - Σ a u.

Parameters:

Name	Type	Description	Default
`X`	`FData`	Curves with the same number of variables and domain as the fitted ones.	required
`forcing`	`FData or sequence`	Forcing functions for `X`, structured as in :meth:`fit`; required exactly when the model was fitted with forcing.	`None`

Returns:

Type	Description
`FData`	The residual functions, in the basis of `X`.

Raises:

Type	Description
`TypeError`	If `X` is not an :class:`~fdatools.core.FData`.
`ValueError`	If `X` or `forcing` does not match the fitted model.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> pda = PDA(order=2, n_grid=None).fit(fd)
>>> new = fdt.FData(np.array([0.0, 2.0, -1.0]), basis)
>>> bool(np.max(np.abs(pda.transform(new)(np.linspace(0, 6, 7)))) < 1e-10)
True

fit_transform ¶

fit_transform(
    X: FData, y: Any = None, *, forcing: Any = None
) -> FData

Fit to X (with its forcing) and return the residual functions.

Parameters:

Name	Type	Description	Default
`X`	`FData`	The curves.	required
`y`	`None`	Ignored; present for scikit-learn compatibility.	`None`
`forcing`	`FData or sequence`	Forcing functions, as in :meth:`fit`.	`None`

Returns:

Type	Description
`FData`	`L x - Σ a u` for the fitted curves (equal to `residuals_`).

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> PDA(order=2).fit_transform(fd).n_curves
2

solve ¶

solve(
    t: Any,
    initial: Any,
    *,
    forcing: Any = None,
    rtol: float = _ODE_RTOL,
    atol: float = _ODE_ATOL,
) -> Array

Integrate the fitted equation from an initial state.

The equation is rewritten as a first-order system in (x, Dx, ..., D^{m-1} x) and integrated with :func:scipy.integrate.solve_ivp (method DOP853) from t[0], forwards or backwards. Without forcing this is the homogeneous equation L x = 0 (the free response, also for a model fitted with forcing); with it, L x = Σ a_k u_k.

Parameters:

Name	Type	Description	Default
`t`	`array_like`	Strictly monotone points inside the domain; the initial condition holds at `t[0]`.	required
`initial`	`array_like`	`(x, Dx, ..., D^{m-1} x)` at `t[0]`: shape `(order,)` for one variable, `(n_vars, order)` for a system.	required
`forcing`	`FData or sequence`	The input functions, structured as in :meth:`fit`, one curve each.	`None`
`rtol`	`float`	Relative and absolute tolerances of the integrator.	`_ODE_RTOL`
`atol`	`float`	Relative and absolute tolerances of the integrator.	`_ODE_RTOL`

Returns:

Type	Description
`ndarray`	`x(t)`: shape `(n_t,)` for one variable, `(n_t, n_vars)` for a system.

Raises:

Type	Description
`ValueError`	If `t` is empty, not monotone or outside the domain, if `initial` has the wrong shape, if `forcing` does not match the fitted model, or if the integrator fails.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> pda = PDA(order=2, n_grid=None).fit(fd)
>>> t = np.linspace(0.0, np.pi, 5)
>>> bool(np.allclose(pda.solve(t, [0.0, 1.0]), np.sin(t), atol=1e-8))
True

stability ¶

stability(
    t: Any = None,
    *,
    n_points: int = _R_GRID,
    forcing: Any = None,
) -> PDAStability

Eigenvalues and equilibria of the fitted equation over time.

Replaces R's eigen.pda. At each time the equations are the first-order system Dz = A(t) z + f(t) in the state z = (x_1, Dx_1, …, D^{m-1} x_1, x_2, …): A(t) has ones on the superdiagonal of each variable's block and -β_{ikj}(t) in the last row of block i, column (k, j); f(t) holds Σ_k a_k(t) u_k(t) in the last row of each block. See :class:PDAStability.

Parameters:

Name	Type	Description	Default
`t`	`array_like`	Time points in the domain. Defaults to `n_points` equally spaced points spanning it (R's `argvals`).	`None`
`n_points`	`int`	Number of default time points. Default `501`, as in R.	`_R_GRID`
`forcing`	`FData or sequence`	The input functions, structured as in :meth:`fit`, one curve each. Defaults to the forcing the model was fitted with when every one of those functions has a single curve.	`None`

Returns:

Type	Description
`PDAStability`	Eigenvalues and equilibrium states at every time point.

Raises:

Type	Description
`ValueError`	If `t` is invalid, or the forcing is missing or does not match the fitted model.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> from fdatools.smoothing import smooth
>>> t = np.linspace(0.0, 1.0, 101)
>>> basis = fdt.BSpline(domain=(0.0, 1.0), n_basis=24, order=5)
>>> fd = smooth(0.5 * (1 - np.exp(-4 * t)), t, basis=basis, lam=0.0).fd
>>> u = fdt.FData(np.array([1.0]), fdt.Constant(domain=(0.0, 1.0)))
>>> result = PDA(order=1).fit(fd, forcing=u).stability(n_points=3)
>>> np.round(result.eigenvalues.real, 6) + 0.0, np.round(result.limits, 6) + 0.0
(array([[-4.],
       [-4.],
       [-4.]]), array([[0.5],
       [0.5],
       [0.5]]))

plot_overlay ¶

plot_overlay(
    ax: Axes | None = None,
    *,
    n_points: int = _PLOT_POINTS,
    labels: Mapping[float, str] | None = None,
    **kwargs: Any,
) -> Axes

Plot the weights of a second-order equation on its stability diagram.

Replaces R's pda.overlay. The trajectory (β₁(t), β₀(t)) is drawn over the parabola β₀ = β₁² / 4: above it the characteristic roots of D²x + β₁ Dx + β₀ x = 0 are complex (the system oscillates), below it they are real; β₁ > 0 damps and β₁ < 0 amplifies.

Parameters:

Name	Type	Description	Default
`ax`	`Axes`	Axes to draw on. A new figure is created when omitted.	`None`
`n_points`	`int`	Number of time points on the trajectory.	`_PLOT_POINTS`
`labels`	`mapping of float to str`	Text to write at the trajectory's position at given times.	`None`
`**kwargs`	`Any`	Passed to :meth:`matplotlib.axes.Axes.plot` for the trajectory.	`{}`

Returns:

Type	Description
`Axes`	The axes drawn on; `lines[0]` is the trajectory and `lines[1]` the parabola.

Raises:

Type	Description
`ValueError`	If the model is not a single second-order equation.

Examples:

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import PDA
>>> basis = fdt.Fourier(domain=(0.0, 2 * np.pi), n_basis=3)
>>> fd = fdt.FData(np.array([[0.0, 0.0], [1.0, 0.3], [0.2, 1.0]]), basis)
>>> len(PDA(order=2).fit(fd).plot_overlay().lines) >= 2
True

phase_plane ¶

phase_plane(
    fd: FData,
    t: Any = None,
    *,
    deriv: tuple[int, int] = (1, 2),
    labels: Mapping[float, str] | None = None,
    ax: Axes | None = None,
    **kwargs: Any,
) -> Axes

Plot one derivative of each curve against another.

Replaces R's phaseplanePlot. The default plots acceleration D²x against velocity Dx: kinetic energy grows with the horizontal distance from the origin and potential energy with the vertical one, so an oscillation traces a loop around the origin.

Parameters:

Name	Type	Description	Default
`fd`	`FData`	Single-variable curves.	required
`t`	`array_like`	Evaluation points. Defaults to 201 equally spaced points on the domain.	`None`
`deriv`	`tuple of int`	Derivative orders on the horizontal and vertical axes. Default `(1, 2)`.	`(1, 2)`
`labels`	`mapping of float to str`	Text written on every curve at its position at the given times, e.g. month names on a year of data.	`None`
`ax`	`Axes`	Axes to draw on. A new figure is created when omitted.	`None`
`**kwargs`	`Any`	Passed to :meth:`matplotlib.axes.Axes.plot`.	`{}`

Returns:

Type	Description
`Axes`	The axes drawn on, with one line per curve.

Raises:

Type	Description
`ValueError`	If `fd` has more than one variable.

Examples:

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.dynamics import phase_plane
>>> basis = fdt.Fourier(domain=(0.0, 1.0), n_basis=5)
>>> fd = fdt.FData(np.array([[0.0, 1.0, 0.0, 0.2, 0.0]]).T, basis)
>>> len(phase_plane(fd).lines)
1

Dynamics¶

PDAStability dataclass ¶

plot ¶

PDA ¶

fit ¶

transform ¶

fit_transform ¶

solve ¶

stability ¶

plot_overlay ¶

phase_plane ¶

PDAStability `dataclass` ¶