Registration

fdatools.registration: continuous registration register(), landmark registration (register(..., landmarks=...) or landmark_register()), the amplitude/phase decomposition RegistrationResult.decompose(), and the scikit-learn Registrator transformer.

Curve registration: remove phase variation by warping the time axis.

A single entry point, :func:register, replaces R fda's register.fd (continuous registration), landmarkreg (landmark registration), register.newfd and AmpPhaseDecomp (via :meth:RegistrationResult.decompose).

Every warping function is strictly increasing by construction. It is built from an unconstrained latent function W exactly as a monotone smooth is (see :func:fdatools.smoothing.smooth with constraint="morph"):

.. math::

h(t) = a + (b - a)\,\frac{\int_a^t e^{W(u)}\,du}{\int_a^b e^{W(u)}\,du},
\qquad W(u) = \sum_k c_k \varphi_k(u),

so h(a) = a, h(b) = b and h' > 0. The registered curve is x(h(t)) (x(h(t) + δ) with a per-curve shift δ for periodic data).

Continuous registration chooses c for each curve by minimising, on the uniform grid t_1 .. t_n with n = max(201, 10 K + 1) points (K the size of the curves' basis), either the least-squares criterion

.. math:: F_1 = \frac1n \sum_j \big(x_0(t_j) - x(h(t_j))\big)^2

or the minimum-eigenvalue criterion of Ramsay & Silverman (2005, §7.6),

.. math::

F_2 = 2\,\mu_{\min}\!\begin{pmatrix}
    \overline{x_0^2} & \overline{x_0\,x\circ h} \\
    \overline{x_0\,x\circ h} & \overline{(x\circ h)^2}
\end{pmatrix},

where the bar is the grid mean, plus the roughness penalty λ cᵀRc. These are the discretisations R fda 6.3.0 uses, so lam means the same thing in both.

Multivariate curves (n_vars > 1, such as the hip and knee angles of the gait data, or the x and y coordinates of handwriting) get one warp per curve, shared by all variables: the fitting criterion is the weighted sum Σ_v w_v F(x_{0v}, x_v∘h) of the criteria of the variables (Ramsay & Silverman 2005, §7.6), with equal weights by default. R fda 6.3.0's register.fd accepts multivariate curves but fits their warps to the first variable alone (measured, black box); var_weights=[1, 0, ...] reproduces that. The first coefficient is pinned to zero, which removes the one exactly flat direction (W -> W + s leaves h unchanged). The minimisation is a safeguarded Newton iteration with the exact Hessian: the first and second derivatives of h with respect to c are cumulative integrals of φ_k e^W and φ_k φ_l e^W, computed by Gauss-Legendre quadrature.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=15)
>>> t = np.linspace(0.0, 1.0, 400)
>>> shifts = [0.42, 0.5, 0.58]
>>> curves = np.stack([np.exp(-((t - s) / 0.1) ** 2) for s in shifts], axis=1)
>>> fd = FData(np.linalg.lstsq(basis(t), curves, rcond=None)[0], basis)
>>> wbasis = BSpline(domain=(0.0, 1.0), n_basis=5)
>>> res = register(fd, criterion="least_squares", warp_basis=wbasis, lam=1e-3)
>>> peaks = t[np.argmax(res.registered(t), axis=0)]
>>> bool(np.ptp(peaks) < 0.5 * np.ptp(shifts))
True

AmpPhaseDecomposition

Bases: NamedTuple

Amplitude/phase split of the variation in a set of curves.

Replaces the list returned by R's AmpPhaseDecomp. Being a named tuple it unpacks as amp_mse, phase_mse, rsq, c = res.decompose().

Attributes:

Name	Type	Description
amp_mse	float	Mean squared error due to amplitude variation (R MS.amp).
phase_mse	float	Mean squared error due to phase variation (R MS.pha). It can be negative when the registration does not improve the alignment.
rsq	float	Share of the total variation that is due to phase, phase_mse / (amp_mse + phase_mse) (R RSQR).
c	float	The constant of the decomposition (R C); it equals one when the warp derivatives are uncorrelated with the squared registered curves.

Examples:

>>> from fdatools.registration import AmpPhaseDecomposition
>>> AmpPhaseDecomposition(1.0, 3.0, 0.75, 1.0).rsq
0.75

RegistrationResult dataclass

RegistrationResult(
    registered: FData,
    warp: FData,
    unregistered: FData,
    latent: FData,
    shift: Array,
    target: FData | None = None,
    warp_inverse: FData | None = None,
    criterion: Array | None = None,
    n_iter: Array | None = None,
)

Outcome of :func:register.

Attributes:

Name	Type	Description
registered	FData	The registered curves x_i(h_i(t)) (R regfd), in the basis of the input curves: the least-squares fit to the warped values on the registration grid.
warp	FData	The warping functions h_i (R warpfd). For continuous registration they are expressed in the warp basis and include the shift of a periodic registration; for landmark registration they are expressed in the basis of the curves. Both are least-squares fits on the registration grid, as in R. :meth:warp_values evaluates the exact warps.
unregistered	FData	The input curves (R yfd).
latent	FData	The latent functions W_i that define the warps (R Wfd).
shift	array	Per-curve time shift δ_i of a periodic registration; zeros otherwise.
target	FData or None	The target curve(s) of a continuous registration; None for landmarks.
warp_inverse	FData or None	The inverse warps h_i⁻¹ of a landmark registration (R warpinvfd), in the basis of the curves; None for continuous registration.
criterion	array or None	Final value of the fitting criterion per curve (continuous only).
n_iter	array or None	Newton iterations used per curve (continuous only).

Notes

For multivariate curves registered, unregistered and target keep their variables (coefficients (n_basis, n_curves, n_vars)), while warp, latent, shift and warp_inverse hold one warp per curve, shared by its variables.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=8)
>>> fd = FData(np.random.default_rng(0).standard_normal((8, 3)), basis)
>>> res = register(fd, landmarks=[[0.4], [0.5], [0.6]])
>>> res.registered.n_curves, res.warp_inverse.n_curves
(3, 3)

warp_values

warp_values(t: Any) -> Array

Evaluate the exact warping functions h_i(t) + δ_i.

Parameters:

Name	Type	Description	Default
t	array	Points in the domain of the curves, shape (n_points,).	required

Returns:

Type	Description
array	Warp values of shape (n_points, n_curves): a tensor (constant, on the device of the result) for a registration of tensors or for tensor t, a NumPy array otherwise.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=8)
>>> fd = FData(np.random.default_rng(1).standard_normal((8, 2)), basis)
>>> res = register(fd, landmarks=[[0.45], [0.55]])
>>> values = res.warp_values(np.linspace(0.0, 1.0, 5))
>>> bool(np.all(np.diff(values, axis=0) > 0))
True

apply

apply(fd: FData) -> FData

Warp new curves with the warps of this registration.

Replaces R's register.newfd: curve i of fd becomes x_i(h_i(t) + δ_i) with the exact warp h_i and shift δ_i estimated here (warped times wrap around the domain when the shifts of a periodic registration are not all zero). The result is the least-squares fit of those values, on the registration grid of fd's basis, in fd's basis -- the same projection :func:register uses, so res.apply(res.unregistered) reproduces res.registered. A typical use is carrying warps estimated on one set of curves over to derived or companion curves (derivatives, other variables).

Parameters:

Name	Type	Description	Default
fd	FData	Curves to warp, univariate or multivariate, with one curve per warp and on the domain of the warps.	required

Returns:

Type	Description
FData	The warped curves in fd's basis, with fd's variables. With :class:torch.Tensor coefficients the result holds a tensor of the same dtype and device, differentiable with respect to fd.coefs (the warps are constants).

Raises:

Type	Description
ValueError	If fd has a different number of curves than there are warps, lies on another domain, or has NaN or infinite coefficients.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=8)
>>> fd = FData(np.random.default_rng(5).standard_normal((8, 3)), basis)
>>> res = register(fd, landmarks=[[0.4], [0.5], [0.6]])
>>> bool(np.allclose(res.apply(fd).coefs, res.registered.coefs))
True
>>> res.apply(fd.derivative()).n_curves
3

decompose

decompose(
    domain: tuple[float, float] | None = None,
) -> AmpPhaseDecomposition

Split the variation of the curves into amplitude and phase parts.

Replaces R's AmpPhaseDecomp (Kneip & Ramsay 2008). With x_i the unregistered curves, y_i the registered ones and h_i the warps,

.. math::

C = 1 + \frac{\frac{1}{N-1}\sum_i \int (Dh_i - \overline{Dh})
              (y_i^2 - \overline{y^2})}{\frac1N \sum_i \int y_i^2},
\quad
\text{amp} = \frac{C}{N}\sum_i \int (y_i - \bar y)^2,
\quad
\text{phase} = C\int \bar y^2 - \int \bar x^2 ,

and rsq = phase / (amp + phase). The integrals use the trapezoidal rule on max(201, 10 K + 1) equally spaced points, K being the size of the basis of the unregistered curves, which is the rule R uses. For multivariate curves the squares are squared Euclidean norms over the variables (every integral is summed over the variables); R's AmpPhaseDecomp rejects multivariate curves.

Parameters:

Name	Type	Description	Default
domain	tuple of float	Sub-interval over which to decompose. Defaults to the whole domain.	None

Returns:

Type	Description
AmpPhaseDecomposition	(amp_mse, phase_mse, rsq, c).

Raises:

Type	Description
ValueError	If fewer than two curves are registered, or domain is not inside the domain of the curves.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=10)
>>> fd = FData(np.random.default_rng(2).standard_normal((10, 4)), basis)
>>> amp, phase, rsq, c = register(fd, landmarks=[[0.4], [0.45], [0.55], [0.6]]).decompose()
>>> bool(abs(rsq - phase / (amp + phase)) < 1e-12)
True

Registrator

Registrator(
    warp_basis: BSpline | None = None,
    *,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    criterion: str = "eigen",
    periodic: bool = False,
    max_iter: int = 100,
    tol: float = 1e-10,
    basis: Basis | None = None,
)

Bases: TransformerMixin, BaseEstimator

Continuous registration as a scikit-learn transformer.

fit registers the training curves to their mean and keeps that mean as the target; transform registers new curves to the same target and returns the coefficients of the registered curves, one row per curve.

Parameters:

Name	Type	Description	Default
warp_basis	BSpline	Basis of the latent warp functions, as in :func:register.	None
lam	float	Roughness weight on the warps. Default 0.	0.0
penalty	int or LDO	Roughness operator on the warps. Default 2.	2
criterion	('eigen', 'least_squares')	Registration criterion. Default "eigen".	"eigen"
periodic	bool	Estimate a periodic time shift as well. Default False.	False
max_iter	int	Newton iterations per curve. Default 100.	100
tol	float	Convergence tolerance. Default 1e-10.	1e-10
basis	Basis	Basis to attach to a plain coefficient matrix; ignored for FData input. None builds a cubic B-spline on (0, 1) with one function per column.	None

Notes

Being a scikit-learn transformer it works on univariate curves (one row of coefficients per curve); multivariate FData raise ValueError.

Attributes:

Name	Type	Description
target_	FData	The registration target (mean of the training curves).
result_	RegistrationResult	The registration of the training curves.
n_features_in_	int	Number of basis coefficients per curve seen during fit.
n_iter_	int	Largest number of Newton iterations used by a training curve.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import Registrator
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=9)
>>> fd = FData(np.random.default_rng(4).standard_normal((9, 5)), basis)
>>> Registrator(criterion="least_squares").fit_transform(fd).shape
(5, 9)

fit

fit(X: Any, y: Any = None) -> Registrator

Register the training curves to their mean.

Parameters:

Name	Type	Description	Default
X	FData or array of shape (n_samples, n_basis)	Curves, or their basis coefficients one curve per row.	required
y	ignored	Present for API compatibility.	None

Returns:

Type	Description
Registrator	The fitted estimator.

transform

transform(X: Any) -> Array

Register curves to the fitted target and return their coefficients.

Parameters:

Name	Type	Description	Default
X	FData or array of shape (n_samples, n_basis)	Curves in the basis seen during fit.	required

Returns:

Type	Description
ndarray	Registered coefficients of shape (n_samples, n_basis).

register

register(
    fd: FData,
    target: FData | None = None,
    *,
    landmarks: Any = None,
    target_landmarks: Any = None,
    warp_basis: Basis | None = None,
    lam: float | None = None,
    penalty: int | LDO = 2,
    criterion: str = "eigen",
    periodic: bool = False,
    init: Any = None,
    init_shift: Any = None,
    max_iter: int = 100,
    tol: float = 1e-10,
    var_weights: Any = None,
) -> RegistrationResult

Register curves by warping their time axis.

Without landmarks this is continuous registration (R register.fd): each curve is warped to minimise its distance to target (the mean curve by default). With landmarks it is landmark registration (R landmarkreg): each curve is warped so that its landmark times move to target_landmarks (the mean landmark times by default).

Multivariate curves get one warp per curve, shared by all their variables; the continuous criterion is then the sum over the variables, weighted by var_weights.

Parameters:

Name	Type	Description	Default
fd	FData	The curves to register, univariate or multivariate.	required
target	FData	One target curve, or one per curve, with the variables of fd. Defaults to fd.mean(). Ignored for landmark registration.	None
landmarks	array	Landmark times, shape (n_curves, n_landmarks) (or (n_curves,) for one landmark), strictly inside the domain and increasing per curve.	None
target_landmarks	array	Where the landmarks should land, shape (n_landmarks,). Defaults to the mean landmark times.	None
warp_basis	Basis	B-spline basis for the latent functions W. Defaults to the linear two-function basis for continuous registration (R's default) and to a cubic spline with knots at the target landmarks otherwise.	None
lam	float	Roughness weight on W. Defaults to 0 for continuous registration (R's default) and to 1e-4 for landmark registration, where the latent coefficients outnumber the landmarks.	None
penalty	int or LDO	Roughness operator on W. Default 2.	2
criterion	('eigen', 'least_squares')	Continuous criterion: the minimum-eigenvalue criterion (R crit=2, the default) or least squares (R crit=1).	"eigen"
periodic	bool	Treat the curves as periodic over their domain and estimate a time shift per curve as well as the warp (R periodic=TRUE).	False
init	FData or array	Starting latent functions: an FData on the warp basis (whose basis then is the warp basis), or a (K,) / (K, n_curves) coefficient array. Defaults to zero, the identity warp.	None
init_shift	float or array	Starting shifts of a periodic registration. Default zero.	None
max_iter	int	Maximum Newton iterations per curve. 0 returns the warps defined by init without optimising. Default 100.	100
tol	float	Convergence tolerance on the largest gradient component, relative to max(1, criterion). Default 1e-10.	1e-10
var_weights	array	Non-negative weight of each variable in the continuous criterion of multivariate curves, shape (n_vars,), at least one positive. Defaults to ones: the criteria of the variables are summed. R fda 6.3.0's register.fd fits multivariate warps to the first variable only, which var_weights=[1, 0, ...] reproduces. Ignored for landmark registration.	None

Returns:

Type	Description
RegistrationResult	The registered curves, warps, latent functions and shifts. With PyTorch input (see Notes) every array and coefficient matrix in it is a tensor of the input's dtype and device.

Raises:

Type	Description
ValueError	On a non-B-spline warp basis, badly shaped landmarks or targets, a target whose number of variables differs from fd's, invalid var_weights, an unknown criterion, a negative lam, or NaN or infinite values in fd, target, init or init_shift.

Warns:

Type

Description

RuntimeWarning

If some curve has not converged after max_iter iterations. This includes a criterion without a finite minimum: with lam=0 the eigenvalue criterion can keep decreasing as the warp squeezes a curve into a vanishing stretch of the domain, and W grows without bound. The line search rejects steps at which the criterion or its derivatives overflow and a failed eigendecomposition of the Hessian falls back to steepest descent, so such curves stop at the last finite iterate (an increasing warp) with this warning rather than an exception. A positive lam keeps W bounded.

Notes

PyTorch input. When fd (or, for continuous registration, an FData target) has :class:torch.Tensor coefficients, the result holds tensors in the dtype and on the device of that input, and registered is differentiable with respect to the input coefficients. Continuous registration then evaluates the criterion in PyTorch and takes its gradient and Hessian by automatic differentiation instead of the analytic formulas; the safeguarded Newton iteration, the line search and the stopping rules are the same, so both paths reach the same optimum. The optimisation always runs in float64 on the CPU (a registration needs tol-level accuracy, and some devices have no float64). Gradients flow through the final evaluation x_i(h_i(t)) with the optimal warps held fixed: the dependence of the optimum on the curves (the implicit gradient) is not propagated. Landmark registration needs no optimisation; its warps depend only on the landmarks, so the gradient of registered with respect to the coefficients is exact. warp, latent, warp_inverse, shift, criterion and n_iter are constant tensors. Multivariate tensor curves work the same way. import fdatools never imports torch; this path imports it on first use.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import register
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=12)
>>> t = np.linspace(0.0, 1.0, 300)
>>> curves = np.stack([np.sin(np.pi * t**p) for p in (0.8, 1.0, 1.25)], axis=1)
>>> fd = FData(np.linalg.lstsq(basis(t), curves, rcond=None)[0], basis)
>>> peaks = t[np.argmax(curves, axis=0)]
>>> res = register(fd, landmarks=peaks)
>>> bool(np.allclose(res.warp_values(np.full(1, peaks.mean())), peaks, atol=1e-3))
True

Multivariate curves (here two variables that share their phase) get one warp per curve:

>>> both = FData(np.stack([fd.coefs, 2.0 * fd.coefs], axis=2), basis)
>>> multi = register(both, criterion="least_squares", warp_basis=BSpline(n_basis=4))
>>> multi.registered.coefs.shape, multi.latent.coefs.shape
((12, 3, 2), (4, 3))

landmark_register

landmark_register(
    fd: FData,
    landmarks: Any,
    target_landmarks: Any = None,
    **kwargs: Any,
) -> RegistrationResult

Landmark registration; shorthand for register(fd, landmarks=...).

Parameters:

Name	Type	Description	Default
fd	FData	The curves to register, univariate or multivariate (one warp per curve, shared by its variables).	required
landmarks	array	Landmark times, shape (n_curves, n_landmarks) or (n_curves,).	required
target_landmarks	array	Target landmark times; the mean landmark times by default.	None
**kwargs	Any	warp_basis, lam and penalty, as in :func:register.	{}

Returns:

Type	Description
RegistrationResult	As returned by :func:register.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.registration import landmark_register
>>> fd = FData(np.random.default_rng(3).standard_normal((6, 2)), BSpline(n_basis=6))
>>> landmark_register(fd, [0.3, 0.5]).latent.n_curves
2