`CustomFilter`#

class CustomFilter(propagator: Propagator | MissionSequence, measurements: MeasurementSpec, settings: FilterSettings, estimator=None, smoother=None, traj_name: str | None = None, traj_dir: str | None = None, overwrite_traj: bool = True, verbose: bool = True)#

Bases: FilterOD

Estimator whose whole algorithm is supplied from a script.

CustomFilter completes the trio started by CustomDynamicModel and CustomMeasurement. Where those replace one model, this replaces the estimator itself: the user writes the estimation algorithm and Scarabaeus provides the scaffolding around it — reference trajectory, measurement timeline, residuals and partials, state-transition matrices, process noise, and the SolutionOD packaging that the rest of the toolchain consumes.

Unlike LKF or SRIF, nothing here assumes a particular estimator structure. The algorithm is free to be

sequential — any Kalman/information/square-root recursion,
batch — accumulate normal equations over all epochs and solve once,
sigma-point (UKF/CKF) — propagate sigma points through the real nonlinear dynamics and measurements with predict() and predict_measurements(),
Gaussian mixture — run any of the above per component and recombine,
polynomial chaos — evaluate the model at the quadrature/collocation nodes and fit the expansion,
sampling (MCMC, particle) — score candidate states by their measurement likelihood.

The estimator is a callable invoked once per fit() iteration. It receives the filter itself — the toolkit below — and returns the histories the solution is built from.

Toolkit. Read the problem:

`timeline`	epochs the solution is reported on
`measurement_epochs`	epochs that carry an observation
`n_states`	length of the estimated state
`initial_covariance_matrix`	\(P_0\)
`initial_deviation`	\(\bar x_0\)
`measurement_block()`	residuals, partials, sigmas at an epoch
`measurement_covariance()`	\(R_k\) from those sigmas
`stm()`	\(\Phi(t, t_0)\)
`stm_between()`	\(\Phi(t_1, t_0)\)
`process_noise_at()`	discrete \(Q_k\)
`safe_inv()`	pseudo-inverse used across the filters

Step outside the linearisation — these run the real nonlinear models and are what make sigma-point, mixture, chaos and sampling methods possible:

`predict()`	propagate an arbitrary initial state
`predict_measurements()`	computed measurements for such a state

Report the answer with record(), or return a dict with the same keys.

Parameters:

propagatorPropagator or MissionSequence: Reference trajectory generator, as for any filter. A sequence must carry the same state on every leg; one that adds or drops states at a node (e.g. finite-burn parameters) raises CustomFilter.SequenceError.
measurementsMeasurementSpec: Measurement models and observations.
settingsFilterSettings: Filter configuration (initial covariance, process noise, output).
estimatorcallable: The estimation algorithm. Invoked as estimator(filter=self) — declare filter (and optionally timeline, n_states, covariance, deviation) in its signature. Must either call record() or return a mapping with "deviation" and "covariance".
smoothercallable, optional: Backward pass, invoked the same way when fit(if_sequential_smooth=True).
traj_name, traj_dir, overwrite_traj, verbose: Forwarded to FilterOD.

Attributes:

initial_covariance_matrix: A priori covariance \(P_0\), without units.
initial_deviation: A priori state deviation \(\bar x_0\), shape (n, 1).
measurement_epochs: Epochs that carry at least one observation.
n_states: Length of the estimated state.
timeline: Epochs the solution is reported on (measurements together with events).

Methods

`measurement_block`(epoch)	Return the measurement block at epoch.
`measurement_covariance`(sigmas, epoch[, ...])	Build the measurement-noise covariance \(R_k\).
`predict`(initial_state[, epochs])	Propagate an arbitrary initial state through the real dynamics.
`predict_measurements`(initial_state[, epochs])	Computed measurements for an arbitrary initial state.
`process_noise_at`(index)	Discrete process-noise covariance \(Q_k\) at a timeline index.
`record`([deviation, covariance, prefit, postfit])	Report the estimator's results.
`safe_inv`(matrix[, min_svd])	Numerically guarded inverse, shared with the built-in filters.
`stm`(epoch[, index, idx_leg])	State-transition matrix from the initial epoch to epoch.
`stm_between`(epoch_from, epoch_to[, idx_leg])	One-step state-transition matrix between two epochs.

See also

scarabaeus.LKF: Reference sequential implementation to model an algorithm on.
scarabaeus.CustomDynamicModel: User-defined force model.
scarabaeus.CustomMeasurement: User-defined measurement model.

Notes

Residuals and partials returned by measurement_block() are computed against the current reference trajectory, and fit() re-linearises about the updated reference between iterations. An estimator that works in absolute states rather than deviations should use predict() and predict_measurements() instead.

Those two run a full propagation (and, for the second, a full measurement pass) per call. A sigma-point filter therefore costs 2n+1 of them per iteration — acceptable offline, but it is the dominant cost, so cache aggressively inside the estimator.

Examples

A sequential Kalman recursion written from scratch:

import numpy as np
import scarabaeus as scb

def sequential(filter):
    P = filter.initial_covariance_matrix
    x = filter.initial_deviation
    deviations, covariances = [], []

    previous = None
    for epoch in filter.timeline:
        if previous is not None:                      # time update
            Phi = filter.stm_between(previous, epoch)
            x, P = Phi @ x, Phi @ P @ Phi.T
        previous = epoch

        block = filter.measurement_block(epoch)
        if block is not None:                         # measurement update
            H = block["partials"][:, : filter.n_states]
            R = filter.measurement_covariance(block["sigma"], epoch)
            S = H @ P @ H.T + R
            K = P @ H.T @ filter.safe_inv(S)
            x = x + K @ (block["residuals"] - H @ x)
            P = (np.eye(filter.n_states) - K @ H) @ P

        deviations.append(x)
        covariances.append(P)

    filter.record(deviation=deviations, covariance=covariances)

od = scb.CustomFilter(prop, measurements, settings, estimator=sequential)
solution, n_iterations, converged = od.fit(max_iterations=3)

A batch least-squares estimator over the same scaffolding:

def batch(filter):
    n = filter.n_states
    information = filter.safe_inv(filter.initial_covariance_matrix)
    normal = np.zeros((n, 1))

    for epoch in filter.measurement_epochs:
        block = filter.measurement_block(epoch)
        H = block["partials"][:, :n] @ filter.stm(epoch)
        W = filter.safe_inv(filter.measurement_covariance(block["sigma"], epoch))
        information += H.T @ W @ H
        normal += H.T @ W @ block["residuals"]

    P0 = filter.safe_inv(information)
    x0 = P0 @ normal

    deviations, covariances = [], []
    for epoch in filter.timeline:
        Phi = filter.stm(epoch)
        deviations.append(Phi @ x0)
        covariances.append(Phi @ P0 @ Phi.T)

    filter.record(deviation=deviations, covariance=covariances)

measurement_block(epoch) → dict | None#

Return the measurement block at epoch.

Parameters:

epochfloat: Epoch, TDB seconds past J2000.

Returns:

dict or None: Keys include residuals, partials, sigma, instruments, measurement_types, dataset_names and t2. None when the epoch carries no observation (an event node, for instance).

measurement_covariance(sigmas, epoch, meas_corr=None, underweighting_factor: float = 1.0) → ndarray#

Build the measurement-noise covariance \(R_k\).

Parameters:

sigmasarray_like: Per-measurement standard deviations from measurement_block().
epochfloat: Epoch the measurements belong to.
meas_corrlist, optional: Off-diagonal correlation terms.
underweighting_factorfloat, optional: Multiplicative inflation of R. Defaults to 1.0.

Returns:

numpy.ndarray: Measurement covariance (m, m).

predict(initial_state, epochs=None)#

Propagate an arbitrary initial state through the real dynamics.

This is the escape hatch from the linearisation: sigma-point, mixture, chaos and sampling estimators need the nonlinear flow, not Phi.

Parameters:

initial_statearray_like: Initial state values, shape (n,), in the reference state’s own order and units.
epochsEpochArray, optional: Epochs to report. Defaults to the propagator’s own span.

Returns:

numpy.ndarray: Propagated states, shape (len(epochs), n).

Notes

Each call runs a full propagation; cache the results inside the estimator rather than calling it per epoch.

predict_measurements(initial_state, epochs=None) → dict#

Computed measurements for an arbitrary initial state.

Propagates initial_state and evaluates every configured measurement model along the resulting trajectory, giving \(h(x)\) without any linearisation.

Parameters:

initial_statearray_like: Initial state values, shape (n,).
epochsEpochArray, optional: Epochs to evaluate at. Defaults to the propagator’s own span.

Returns:

dict: Mapping of dataset name to the computed measurements.

Notes

This is the most expensive call in the toolkit: one propagation plus a full measurement pass. A sigma-point filter needs 2n+1 per iteration.

process_noise_at(index: int) → ndarray | None#

Discrete process-noise covariance \(Q_k\) at a timeline index.

Parameters:

indexint: Position in timeline.

Returns:

numpy.ndarray or None: Q_k, or None when no process noise is configured.

record(deviation=None, covariance=None, prefit=None, postfit=None)#

Report the estimator’s results.

Parameters:

deviationsequence, optional: State deviation per timeline epoch, each (n, 1).
covariancesequence, optional: Covariance per timeline epoch, each (n, n).
prefit, postfitlist, optional: Residual records in the layout the built-in filters use: [residuals, sigmas, instruments, types, dataset_names, epochs].

safe_inv(matrix, min_svd: float = 1e-14) → ndarray#

Numerically guarded inverse, shared with the built-in filters.

Parameters:

matrixnumpy.ndarray: Matrix to invert.
min_svdfloat, optional: Singular-value floor. Defaults to 1e-14.

Returns:

numpy.ndarray: The (pseudo-)inverse.

stm(epoch, index=None, idx_leg=None) → ndarray#

State-transition matrix from the initial epoch to epoch.

Parameters:

epochfloat: Target epoch, TDB seconds.
indexint, optional: Position of epoch in the timeline. Supplying it lets the trajectory pick the stored matrix directly instead of searching.
idx_legint, optional: Leg index, for sequence trajectories.

Returns:

numpy.ndarray: Phi(epoch, t0), shape (n, n).

stm_between(epoch_from, epoch_to, idx_leg=None) → ndarray#

One-step state-transition matrix between two epochs.

\[\Phi(t_1, t_0) = \Phi(t_1, t_{\text{ref}})\, \Phi(t_0, t_{\text{ref}})^{-1}\]

Parameters:

epoch_from, epoch_tofloat: Bounding epochs, TDB seconds.
idx_legint, optional: Leg index, for sequence trajectories.

Returns:

numpy.ndarray: Phi(epoch_to, epoch_from), shape (n, n).

property initial_covariance_matrix: ndarray#: A priori covariance \(P_0\), without units.

property initial_deviation: ndarray#: A priori state deviation \(\bar x_0\), shape (n, 1).

property measurement_epochs: ndarray#: Epochs that carry at least one observation.

property n_states: int#: Length of the estimated state.

property timeline: ndarray#: Epochs the solution is reported on (measurements together with events).

exception ResultError#

Bases: ValueError

Raised when the estimator does not report a usable result.

exception SequenceError#

Bases: ValueError

Raised when a mission sequence changes the state from one leg to the next.

exception ShapeError#

Bases: ValueError

Raised when a reported history has an unexpected shape.

CustomFilter#

`CustomFilter`#