UtilsOD

class UtilsOD

Bases: object

Utility methods for Orbit Determination solution analysis and manipulation.

All methods are static and operate on raw NumPy arrays or SolutionOD objects. No instantiation is required:

d = scb.UtilsOD.closeness(c1, P1, c2, P2)['d']
c_m, P_m = scb.UtilsOD.merge(c1, P1, c2, P2)

See also

scarabaeus.SolutionOD

OD solution container.

scarabaeus.Utils

General-purpose Scarabaeus utilities.

References

[1]

McElrath, T. P., Watkins, M. M., Portock, B. M., Graat, E. J., Baird, D. T., Wawrzyniak, G. G., Guinn, J. R., Antreasian, P. G., Attiyah, A. A., Baalke, R. C., & Taber, W. L. (2004). Mars Exploration Rovers Orbit Determination Filter Strategy. In AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Providence, Rhode Island, 16–19 August 2004. AIAA Paper 2004-5082. https://doi.org/10.2514/6.2004-5082

Methods

closeness(c1, P1, c2, P2[, alpha_range, ...])	Compute the sigma-ellipsoid closeness measure between two OD solutions.
closeness_matrix(solutions[, labels, verbose])	Compute the N×N pairwise closeness matrix for a list of OD solutions.
mahalanobis(x, c, P)	Compute the Mahalanobis distance between a point and a Gaussian.
merge(c1, P1, c2, P2)	Combine two independent Gaussian estimates via the information filter.
merge_many(solutions)	Combine N independent Gaussian estimates sequentially.
nees(x_true, c_est, P_est)	Compute the Normalized Estimation Error Squared (NEES).
plot_closeness_matrix(D[, labels, scheme, ...])	Visualise an N×N pairwise closeness matrix.
sigma_overlap(c1, P1, c2, P2[, n_sigma])	Check whether two sigma ellipsoids overlap.

static closeness(c1: ndarray, P1: ndarray, c2: ndarray, P2: ndarray, alpha_range: tuple = (1e-06, 1000000.0), n_grid: int = 10000, tol: float = 1e-10) → dict

Compute the sigma-ellipsoid closeness measure between two OD solutions.

Given two Gaussian solutions S1 = (c1, P1) and S2 = (c2, P2), returns the smallest t > 0 such that the t-sigma uncertainty ellipsoids are tangent (AIAA 2004-4982 [1]). A value d < 1 means the solutions are mutually consistent within their 1-sigma bounds.

Parameters:

• c1 (ndarray, shape (n,)) – Mean state vector of solution 1.

• P1 (ndarray, shape (n, n)) – Covariance matrix of solution 1 (symmetric positive definite).

• c2 (ndarray, shape (n,)) – Mean state vector of solution 2.

• P2 (ndarray, shape (n, n)) – Covariance matrix of solution 2 (symmetric positive definite).

• alpha_range ((float, float), optional) – Search range (alpha_min, alpha_max) for the Lagrange multiplier. Defaults to (1e-6, 1e6).

• n_grid (int, optional) – Number of logarithmically spaced grid points for the initial root bracket search. Defaults to 10_000.

• tol (float, optional) – Bisection convergence tolerance for alpha root refinement. Defaults to 1e-10.

Returns:

result – Dictionary with the following keys:

• 'd' – float, closeness in units of sigma.

• 't_sq' – float, d ** 2.

• 'alpha' – float, Lagrange multiplier at the tangency point.

• 'x_tang' – ndarray shape (n,), tangency point in state space.

• 'all_roots' – list of all positive alpha roots found.

• 'all_t_sq' – list of t² values for each root.

Return type:

dict

Notes

The algorithm evaluates the stationarity condition

\[f(\alpha) = \Delta c^\top (P_2 + \alpha P_1)^{-1} (\alpha^2 P_1 - P_2) (P_2 + \alpha P_1)^{-1} \Delta c = 0\]

on a logarithmic grid, brackets sign changes, and refines them with bisection. If no root is found a UserWarning is issued and 'd' is returned as nan.

See also

UtilsOD.closeness_matrix

Pairwise matrix for N solutions.

References

[1]

Woffinden & Geller (2004). AIAA 2004-4982.

Examples

import numpy as np
import scarabaeus as scb

c1 = np.zeros(6)
P1 = np.eye(6) * 4.0
c2 = np.array([1.0, 0.5, 0.2, 0.0, 0.0, 0.0])
P2 = np.eye(6) * 2.0

res = scb.UtilsOD.closeness(c1, P1, c2, P2)
print(f"d = {res['d']:.4f} sigma")

static closeness_matrix(solutions: list, labels: list = None, verbose: bool = True, **closeness_kwargs) → tuple

Compute the N×N pairwise closeness matrix for a list of OD solutions.

Parameters:

• solutions (list of (c, P) or SolutionOD) – Each entry is either a (mean_vector, covariance_matrix) tuple or a SolutionOD object (final epoch is used).

• labels (list of str, optional) – Labels for each solution. Defaults to ["S1", "S2", …].

• verbose (bool, optional) – Print progress for each pair. Defaults to True.

• **closeness_kwargs – Forwarded to closeness() (e.g. alpha_range, n_grid).

Returns:

• D (ndarray, shape (N, N)) – Symmetric matrix where D[i, j] = d(Si, Sj) in sigma units. Diagonal entries are zero.

• labels (list of str) – Labels used (auto-generated when not supplied).

See also

UtilsOD.plot_closeness_matrix

Visualise the returned matrix.

static mahalanobis(x: ndarray, c: ndarray, P: ndarray) → float

Compute the Mahalanobis distance between a point and a Gaussian.

\[d_M = \sqrt{(x - c)^\top P^{-1} (x - c)}\]

Parameters:

• x (ndarray, shape (n,)) – Query point.

• c (ndarray, shape (n,)) – Distribution mean.

• P (ndarray, shape (n, n)) – Covariance matrix (symmetric positive definite).

Returns:

d – Mahalanobis distance in units of sigma.

Return type:

float

static merge(c1: ndarray, P1: ndarray, c2: ndarray, P2: ndarray) → tuple

Combine two independent Gaussian estimates via the information filter.

Given two unbiased estimates of the same quantity with covariances P1 and P2, the optimal linear combination is:

\[P_m = (P_1^{-1} + P_2^{-1})^{-1}, \quad c_m = P_m (P_1^{-1} c_1 + P_2^{-1} c_2)\]

This is equivalent to a Kalman update where one estimate treats the other as a measurement [2]_.

Parameters:

• c1 (ndarray, shape (n,)) – Mean of estimate 1.

• P1 (ndarray, shape (n, n)) – Covariance of estimate 1.

• c2 (ndarray, shape (n,)) – Mean of estimate 2.

• P2 (ndarray, shape (n, n)) – Covariance of estimate 2.

Returns:

• c_merged (ndarray, shape (n,)) – Combined mean.

• P_merged (ndarray, shape (n, n)) – Combined covariance (always ≤ min(P1, P2) element-wise on the diagonal).

Raises:

np.linalg.LinAlgError – If either covariance matrix is singular.

Examples

c_m, P_m = scb.UtilsOD.merge(c1, P1, c2, P2)

static merge_many(solutions: list) → tuple

Combine N independent Gaussian estimates sequentially.

Applies merge() iteratively to the list of solutions, accumulating the combined mean and covariance.

Parameters:

solutions (list of (c, P) or SolutionOD) – Ordered list of solutions. Each entry is a (mean, covariance) tuple or a SolutionOD object (final epoch).

Returns:

• c_merged (ndarray, shape (n,)) – Combined mean.

• P_merged (ndarray, shape (n, n)) – Combined covariance.

Raises:

ValueError – If fewer than two solutions are provided.

Examples

c_m, P_m = scb.UtilsOD.merge_many([sol1, sol2, sol3])

static nees(x_true: ndarray, c_est: ndarray, P_est: ndarray) → float

Compute the Normalized Estimation Error Squared (NEES).

A standard filter consistency metric [2]_:

\[\epsilon = (x_{true} - c_{est})^\top P_{est}^{-1} (x_{true} - c_{est})\]

For a consistent filter, the expected value of NEES equals the state dimension n.

Parameters:

• x_true (ndarray, shape (n,)) – True state vector.

• c_est (ndarray, shape (n,)) – Filter estimated mean.

• P_est (ndarray, shape (n, n)) – Filter estimated covariance.

Returns:

epsilon – NEES value. Expected value for a consistent filter: n.

Return type:

float

static plot_closeness_matrix(D: ndarray, labels: list = None, scheme: str = 'merb', title: str = 'OD Solution Consistency Matrix', show_values: bool = True, figsize: tuple = None, save_path: str = None) → None

Visualise an N×N pairwise closeness matrix.

Parameters:

• D (ndarray, shape (N, N)) – Closeness matrix returned by closeness_matrix().

• labels (list of str, optional) – Axis tick labels. Defaults to ["S1", "S2", …].

• scheme ({'merb', 'mera', 'sigma', 'continuous'}, optional) –

  Colormap scheme:

  • 'merb' (default) — five-band categorical map.

  • 'mera' — two-band binary map at the 0.8 σ boundary.

  • 'sigma' — four-band map at 1, 2, 3 σ boundaries.

  • 'continuous' — RdYlGn_r continuous colormap.

• title (str, optional) – Figure title. Defaults to 'OD Solution Consistency Matrix'.

• show_values (bool, optional) – Annotate upper-triangle cells with numeric values. Defaults to True.

• figsize ((float, float), optional) – Figure size in inches. Auto-scaled from N when None.

• save_path (str, optional) – Save the figure to this path at 300 dpi instead of displaying it.

Return type:

None

static sigma_overlap(c1: ndarray, P1: ndarray, c2: ndarray, P2: ndarray, n_sigma: float = 1.0) → bool

Check whether two sigma ellipsoids overlap.

Two ellipsoids overlap if the sum of their Mahalanobis radii (each evaluated at the midpoint between the means) exceeds the distance between the means [2]_. This is a fast, conservative heuristic; use closeness() for the exact tangency value.

Parameters:

• c1 (ndarray, shape (n,)) – Mean of distribution 1.

• P1 (ndarray, shape (n, n)) – Covariance of distribution 1.

• c2 (ndarray, shape (n,)) – Mean of distribution 2.

• P2 (ndarray, shape (n, n)) – Covariance of distribution 2.

• n_sigma (float, optional) – Number of sigma defining the ellipsoids. Defaults to 1.0.

Returns:

overlap – True if the n_sigma-ellipsoids overlap.

Return type:

bool