Add generic pairwise covariance association features

docs/api-overview.md

@@ -144,13 +144,15 @@ Common starting points include:
 ### `pyrecest.utils`
 
 Reusable utility helpers for assignment, association models, history recording,
-multi-session assignment, generic track-matrix evaluation, and point-set
-registration.
+multi-session assignment, generic track-matrix evaluation, Gaussian/covariance
+pairwise features, and point-set registration.
 
 Common starting points include:
 
 - `murty_k_best_assignments`;
 - `LogisticPairwiseAssociationModel`;
+- `pairwise_mahalanobis_distances`;
+- `pairwise_covariance_shape_components`;
 - `HistoryRecorder`;
 - `solve_multisession_assignment`;
 - `normalize_track_matrix`;

pyproject.toml

@@ -82,6 +82,8 @@ disable = [
   "too-many-branches",
   "too-many-statements",
   "too-many-arguments",
+  "too-many-positional-arguments",
+  "too-many-instance-attributes",
   "duplicate-code",
 ]
 

src/pyrecest/distributions/conditional/sd_half_cond_sd_half_grid_distribution.py

@@ -57,9 +57,11 @@ def __init__(self, grid, grid_values, enforce_pdf_nonnegative=True):
 
     def _check_normalization(self, tol=0.01):
         """Warn if any column is not normalized to 1 over the hemisphere."""
+        # Match HyperhemisphericalGridDistribution.get_manifold_size(), which
+        # uses the grid embedding dimension for the existing quadrature rule.
         hemisphere_surface = 0.5 * (
             AbstractHypersphereSubsetDistribution.compute_unit_hypersphere_surface(
-                self.grid.shape[1] - 1
+                self.grid.shape[1]
             )
         )
         ints = mean(self.grid_values, axis=0) * hemisphere_surface

src/pyrecest/filters/_linear_gaussian.py

@@ -1,4 +1,4 @@
-# pylint: disable=no-name-in-module,no-member
+# pylint: disable=no-name-in-module,no-member,redefined-outer-name
 """Backend-native linear-Gaussian predict/update primitives."""
 
 from pyrecest.backend import (
@@ -45,8 +45,7 @@ def normalized_innovation_squared(innovation, innovation_covariance):
     innovation_dim = innovation.shape[0]
     if innovation_covariance.shape != (innovation_dim, innovation_dim):
         raise ValueError(
-            "innovation_covariance must have shape "
-            "(innovation_dim, innovation_dim)"
+            "innovation_covariance must have shape (innovation_dim, innovation_dim)"
         )
     return transpose(innovation) @ linalg.solve(innovation_covariance, innovation)
 
@@ -261,8 +260,8 @@ def linear_gaussian_update(
     identity = eye(state_dim)
     correction = identity - kalman_gain @ measurement_matrix
     updated_covariance = correction @ covariance @ transpose(correction)
-    updated_covariance = updated_covariance + kalman_gain @ scaled_meas_noise @ transpose(
-        kalman_gain
+    updated_covariance = (
+        updated_covariance + kalman_gain @ scaled_meas_noise @ transpose(kalman_gain)
     )
     updated_covariance = 0.5 * (updated_covariance + transpose(updated_covariance))
 

src/pyrecest/filters/so3_grid_transition.py

@@ -9,8 +9,9 @@
     clip,
     exp,
     linalg,
-    mean,
     ndim,
+    reshape,
+    sum,
     transpose,
 )
 from pyrecest.distributions._so3_helpers import (
@@ -67,8 +68,9 @@ def so3_right_multiplication_grid_transition(
 
     ``exp(kappa * |<q_next, q_current * delta_q>|**2)``.
 
-    The columns are normalized by the S3 upper-hemisphere grid quadrature rule,
-    so ``mean(grid_values[:, j]) * surface(S3+) == 1`` for every column.
+    The columns are normalized by the hyperhemispherical grid quadrature rule
+    used by :class:`HyperhemisphericalGridFilter`, so
+    ``mean(grid_values[:, j]) * manifold_size == 1`` for every column.
     """
 
     if kappa <= 0.0:
@@ -83,15 +85,17 @@ def so3_right_multiplication_grid_transition(
     # Subtracting the per-column maximum keeps the normalization stable for
     # large kappa without changing the normalized conditional density.
     exponents = kappa * inner_products**2
-    scores = exp(exponents - amax(exponents, axis=0, keepdims=True))
+    column_maxima = reshape(amax(exponents, axis=0), (1, exponents.shape[1]))
+    scores = exp(exponents - column_maxima)
 
     manifold_size = (
         0.5
         * AbstractHypersphereSubsetDistribution.compute_unit_hypersphere_surface(
-            quaternion_grid.shape[1] - 1
+            quaternion_grid.shape[1]
         )
     )
-    column_integrals = mean(scores, axis=0, keepdims=True) * manifold_size
+    column_integrals = sum(scores, axis=0, keepdims=True) / scores.shape[0]
+    column_integrals = column_integrals * manifold_size
     density_values = scores / column_integrals
 
     return SdHalfCondSdHalfGridDistribution(
@@ -121,7 +125,9 @@ def _as_quaternion_grid(grid):
 
     quaternion_grid = array(grid, dtype=float)
     if ndim(quaternion_grid) != 2 or quaternion_grid.shape[1] != 4:
-        raise ValueError("grid must have shape (n_grid, 4) with scalar-last quaternions.")
+        raise ValueError(
+            "grid must have shape (n_grid, 4) with scalar-last quaternions."
+        )
     if quaternion_grid.shape[0] == 0:
         raise ValueError("grid must contain at least one quaternion.")
     if not all(linalg.norm(quaternion_grid, axis=1) > 0.0):

src/pyrecest/filters/so3_product_particle_filter.py

@@ -1,5 +1,6 @@
 """Particle filter for Cartesian products of SO(3)."""
 
+import math
 from collections.abc import Callable
 
 # pylint: disable=no-name-in-module,no-member
@@ -157,7 +158,7 @@ def _normalize_log_weights(log_weights):
     @staticmethod
     def _validate_probability(name: str, value: float) -> float:
         value = float(value)
-        if not isfinite(value) or value < 0.0 or value > 1.0:
+        if not math.isfinite(value) or value < 0.0 or value > 1.0:
             raise ValueError(f"{name} must be a finite probability in [0, 1].")
         return value
 
@@ -198,7 +199,9 @@ def _as_component_noise_std(noise_std, component_noise_std, num_rotations):
 
         if ndim(sigma) == 0:
             if not isfinite(sigma) or sigma <= 0.0:
-                raise ValueError("noise standard deviations must be positive and finite.")
+                raise ValueError(
+                    "noise standard deviations must be positive and finite."
+                )
             return ones(num_rotations) * sigma
         if sigma.shape != (num_rotations,):
             raise ValueError("component_noise_std must have shape (num_rotations,).")
@@ -235,13 +238,15 @@ def confidence_to_noise_std(
         min_sigma = float(noise_std)
         max_sigma = float(max_noise_std)
         exponent = float(confidence_exponent)
-        if min_sigma <= 0.0 or not isfinite(min_sigma):
+        if min_sigma <= 0.0 or not math.isfinite(min_sigma):
             raise ValueError("noise_std must be positive and finite.")
-        if max_sigma <= 0.0 or not isfinite(max_sigma):
+        if max_sigma <= 0.0 or not math.isfinite(max_sigma):
             raise ValueError("max_noise_std must be positive and finite.")
         if max_sigma < min_sigma:
-            raise ValueError("max_noise_std must be greater than or equal to noise_std.")
-        if exponent <= 0.0 or not isfinite(exponent):
+            raise ValueError(
+                "max_noise_std must be greater than or equal to noise_std."
+            )
+        if exponent <= 0.0 or not math.isfinite(exponent):
             raise ValueError("confidence_exponent must be positive and finite.")
 
         variance = min_sigma * min_sigma + (1.0 - confidence) ** exponent * (

src/pyrecest/utils/__init__.py

@@ -22,6 +22,10 @@
     estimate_thin_plate_spline,
     joint_tps_registration_assignment,
 )
+from .pairwise_covariance_features import (
+    pairwise_covariance_shape_components,
+    pairwise_mahalanobis_distances,
+)
 from .track_evaluation import (
     complete_track_set,
     normalize_track_matrix,
@@ -56,6 +60,8 @@
     "estimate_thin_plate_spline",
     "joint_tps_registration_assignment",
     "murty_k_best_assignments",
+    "pairwise_covariance_shape_components",
+    "pairwise_mahalanobis_distances",
     "complete_track_set",
     "normalize_track_matrix",
     "pairwise_track_set",