Metadata-Version: 2.5
Name: gri-convolve
Version: 0.6.1
Summary: Ellipsoid convolution for combining geolocation estimates with outlier detection
Project-URL: Homepage, https://geosolresearch.com
Project-URL: Repository, https://gitlab.com/geosol-foss/python/gri-convolve
Project-URL: Issues, https://gitlab.com/geosol-foss/python/gri-convolve/-/issues
Project-URL: Changelog, https://gitlab.com/geosol-foss/python/gri-convolve/-/releases
Author-email: GeoSol Research Inc <contact@geosolresearch.com>
License-Expression: MIT
License-File: LICENSE
Keywords: clustering,convolution,ellipsoid,fusion,geolocation
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: GIS
Requires-Python: >=3.12
Requires-Dist: gri-ell>=0.2.0
Requires-Dist: gri-pos>=0.2.0
Requires-Dist: gri-utils>=0.5.1
Requires-Dist: numpy>=2.3.3
Requires-Dist: scipy>=1.16.3
Description-Content-Type: text/markdown

[![GeoSol Research Logo](https://geosolresearch.com/logos/foss_logo.png "GeoSol Research")](https://geosolresearch.com)

# Convolve (Ellipsoid Fusion)

Ellipsoid convolution functions for combining geolocation estimates with outlier detection and multi-cluster support.

## Overview

gri-convolve provides three functions of increasing sophistication for fusing collections of `Ell` (ellipsoid) objects into combined position estimates:

- **`convolve`** -- combine all input ellipsoids into a single fused result with no outlier rejection
- **`smart_convolve`** -- iteratively remove outliers by Mahalanobis distance before fusing
- **`cluster_convolve`** -- find multiple clusters within a dataset and fuse each independently

Each function operates on `Ell` objects from gri-ell, which pair a 3D position with a statistical covariance (or information matrix). The output is one or more fused `Ell` objects representing the combined position estimate and its uncertainty.

Recursive single-target tracking -- the `IMM` / `SmartSegmentedIMM` filters, the
motion-model bank, EKF/UKF observable updates, and RTS smoothing -- is provided
by the companion **`gri-kalman`** package. Convolution is the batch face and the
filter is the recursive face of the same estimation problem: the convolver seeds
and refines tracks, the filter does the online sequential update. See the
`gri-kalman` README for the tracker API and EKF-vs-UKF guidance.

Requires Python 3.12+.

## Mathematical Background

**Information matrix fusion.** Given N ellipsoids, each with position `x_k` and information matrix `I_k` (the inverse of the covariance matrix, in XYZ coordinates, 1/m^2, 1-sigma), the fused position and information matrix are:

    S = sum(I_k)              (combined information matrix)
    x = S^{-1} sum(I_k x_k)  (fused position)

This is the maximum-likelihood estimator under Gaussian assumptions.

**Inflation methods.** The raw fusion above underestimates uncertainty when inputs are inconsistent. Three modes control how the output covariance is inflated:

- `"none"` -- strict information matrix combination (no inflation); calibrated when the inputs' own covariances are
- `"std"` -- scale by the inputs' Mahalanobis scatter and add the sample covariance of their offsets
- `"bart"` (default) -- as `"std"`, but each input's horizontal offset is projected onto that input's own minor axis first, so scatter along a major axis, where it is expected, adds nothing

Both inflating modes are conservative on consistent inputs (about 3.5x and 4.6x the calibrated variance for `"bart"` and `"std"` in a Monte Carlo; see the shipped `overview.md`).

**Outlier detection.** `smart_convolve` measures the fused point in each input's own metric:

    d_k = sqrt((x_k - mu)^T I_k (x_k - mu) / 7.815)

where `mu` is the (uninflated) fused position, `I_k` is input k's own information matrix, and 7.815 is the 3-degree-of-freedom 95% chi-square, so `d_k = 1` is the surface of input k's 95% ellipsoid. Inputs with `d_k > max_norm` are removed worst first, refitting each time.

Reference: Mahalanobis, P.C. (1936). "On the generalized distance in statistics."

## Documentation

The wheel ships its documentation inside the package, in `gri_convolve/docs/`,
so it is available wherever the package is installed:

- `overview.md` -- the three functions, the fusion and its inflation modes,
  rejection, clustering, the array fast path, and input checks
- `api_summary.md` -- every public function and signature (generated)

Print the directory with
`python -c "import gri_convolve, pathlib; print(pathlib.Path(gri_convolve.__file__).parent / 'docs')"`.
Every example in those files is run by the test suite.

## Installation

```bash
pip install gri-convolve
```

For development:

```bash
git clone https://gitlab.com/geosol-foss/python/gri-convolve.git
cd gri-convolve
uv sync
```

## Quick Start

```python
from gri_convolve import convolve, smart_convolve, cluster_convolve
from gri_ell import Ell
from gri_pos import Pos
import numpy as np

# Create some ellipsoids at nearby positions (an altitude bound is required)
e1 = Ell.from_2d(Pos.LLA(40.0, -105.0, 1600), 100, 50, 45, alt_95_m=200)
e2 = Ell.from_2d(Pos.LLA(40.001, -105.001, 1610), 120, 60, 30, alt_95_m=200)
e3 = Ell.from_2d(Pos.LLA(40.0005, -104.999, 1605), 90, 45, 50, alt_95_m=200)

# Simple fusion
fused = convolve([e1, e2, e3])
print(fused.lla)            # Fused position
print(fused.ellipse.sma_95) # Fused semi-major axis (95%, meters)
```

## `convolve()`

Fuses all input ellipsoids into a single result. No outlier detection.

```python
fused = convolve(ells, inflation="bart")
```

**Parameters:**

- `ells` -- sequence or generator of `Ell` objects
- `inflation` -- `"none"`, `"std"`, or `"bart"` (default: `"bart"`)

**Returns:** A single fused `Ell`.

## `smart_convolve()`

Fuses with iterative outlier rejection. Computes the fused point, finds the input with the largest normalized Mahalanobis distance, and removes it if it exceeds `max_norm`. Repeats until all remaining inputs are within tolerance or fewer than `min_pts` remain.

```python
result = smart_convolve(ells, max_norm=2.0, min_pts=3)
if result is not None:
    fused_ell, used_indices, discarded_indices = result
```

**Parameters:**

- `ells` -- sequence or generator of `Ell` objects
- `max_norm` -- maximum allowed normalized distance (default: 2.0)
- `min_pts` -- minimum inputs required for a valid result (default: 3)
- `check_scale` -- warn on a likely unit mismatch or near-total rejection (default: True; see [Unit sanity checks](#unit-sanity-checks))

**Returns:** `(Ell, list[int], list[int])` or `None` if no valid cluster is found.

Pre-cluster your data before calling `smart_convolve`. Without pre-clustering, a large group of scattered noise points can cause valid clusters to be discarded first.

## `cluster_convolve()`

Finds multiple clusters within a dataset by peeling them off densest first (`batch=True`, the default) or one at a time with `smart_convolve` (`batch=False`, the exact reference path). Leftover inputs go to the discard list; the used lists and the discard list together hold every input index exactly once.

```python
locations, used_per_location, discarded = cluster_convolve(
    ells,
    max_norm=2.0,
    min_pts=3,
    max_pts=10,
    min_sma_m=50.0,
)

for loc, indices in zip(locations, used_per_location):
    print(f"Cluster at {loc.lla} using {len(indices)} inputs")
```

**Parameters:**

- `ells` -- sequence or generator of `Ell` objects
- `max_norm` -- maximum normalized Mahalanobis distance (default: 2.0)
- `min_pts` -- minimum inputs per cluster (default: 3)
- `max_pts` -- maximum inputs per cluster; splits larger groups (default: None)
- `min_sma_m` -- minimum semi-major axis for output ellipsoids in meters (default: 0)
- `max_ori_spread` -- sort by orientation before splitting for diversity (default: True)
- `check_scale` -- warn on a likely unit mismatch (default: True; see [Unit sanity checks](#unit-sanity-checks))
- `min_cluster_size`, `min_cluster_frac`, `min_utilization` -- size floor for returned clusters, and a coverage floor that keeps peeling (batch path only)
- `alt_post_process` -- callback for altitude correction (e.g., snap to terrain); `gri_convolve.altitude.nearest` is provided

**Returns:** `(list[Ell], list[list[int]], list[int])`

## Unit sanity checks

The convolution math is scale-invariant -- a Mahalanobis distance is dimensionless, so entering every quantity in kilometers instead of meters changes nothing. Trouble appears only when the *position* scale and the *covariance* scale disagree. The common mistake is positions in meters (positions always resolve to ECEF meters) paired with covariance sigmas typed in kilometers: the information matrices come out about a million times too large, every residual distance is about a thousand times too large, so `smart_convolve` and `cluster_convolve` reject nearly every point one at a time -- slow, and with an almost-empty result.

To catch this, `smart_convolve` and `cluster_convolve` run two cheap checks on the full input (controlled by `check_scale`, default `True`) and emit a `ConvolveScaleWarning` when either trips:

- **Scale pre-flight** (before the rejection loop): if the tightest-packed inputs still sit many 1-sigma widths from their nearest neighbor -- meaning no cluster can form -- it warns that the sigmas look too small for the point spacing, and flags a ratio near 1000 as a likely kilometer/meter mix-up. Because it measures nearest-neighbor spacing (a local quantity), legitimate multi-site data with many separate clusters does not trip it.
- **Utilization backstop** (after): if nearly all inputs land in the discard list, it warns that most of the data was rejected. (`cluster_convolve` suppresses this when a cluster size floor is set, since trimming small clusters discards points by design.)

The checks never change the result -- they only warn. Filter or catch them with the exported `ConvolveScaleWarning` category, or pass `check_scale=False` to silence them (e.g. for data that is intentionally spread far in sigma units):

```python
import warnings
from gri_convolve import ConvolveScaleWarning, smart_convolve

with warnings.catch_warnings():
    warnings.simplefilter("error", ConvolveScaleWarning)  # promote to an exception
    result = smart_convolve(ells)
```

## Degenerate inputs are rejected

Separate from the warnings above, and not optional. Every semi-axis of an input ellipsoid must be finite and greater than zero. A zero-width axis -- most often `alt_95_m=0`, written when only a 2D ellipse is available -- makes the covariance rank deficient, and there is no fused answer to give.

This used to fail in whichever of four ways the arithmetic happened to land in, decided by float noise in the ENU-to-ECEF rotation at that latitude and longitude: a bare `LinAlgError: Singular matrix` at build time, the same error much later from the inflation step, a silent `None` with every point discarded, or a silently wrong answer built from a subset. It now raises `ConvolveInputError` up front, naming the offending rows:

```python
from gri_convolve import ConvolveInputError, smart_convolve

try:
    result = smart_convolve(ells)
except ConvolveInputError as exc:
    print(exc)
    # 1 of 5000 input covariances have a variance of zero or less on the
    # diagonal: rows [4999]. A zero-width axis (commonly alt_95_m=0) has no
    # invertible covariance; the convolve solves in 3D, so even a pinned
    # altitude needs a real uncertainty.
```

`ConvolveInputError` subclasses `numpy.linalg.LinAlgError` (itself a `ValueError`), so code already guarding a convolve call against a singular matrix keeps working and simply gets a message that says which input is at fault.

No default is substituted for a zero altitude uncertainty. The solver works in 3D, and whether the right value is a meter, ten meters, or a fraction of the semi-major axis is a modeling decision that belongs to the caller. A 2D-only `Ell` built without `alt_95_m` is rejected with its own message for the same reason.

Eccentricity is not degeneracy: a semi-major/semi-minor ratio of 1e7 is accepted. `params_to_arrays` builds the information matrix directly from the ellipse parameters instead of inverting a covariance, so it keeps full float64 precision at any ratio, where inverting loses accuracy as the square of the ratio.

## Units and Conventions

- Positions are in ECEF XYZ (meters) internally
- Information matrices are in XYZ, 1/m^2, 1-sigma
- Covariance matrices are in ENU, m^2, 1-sigma
- Output ellipse parameters (SMA, SMI, orientation) are at 95% confidence
- Mahalanobis distances are normalized to 95% scale for `max_norm` comparisons

## Dependencies

- **gri-ell**: Ellipsoid objects with position and covariance
- **gri-pos**: Position objects (XYZ, LLA coordinates)
- **gri-utils**: Coordinate conversions and constants
- **numpy**: Array operations
- **scipy**: Nearest-neighbor search (`KDTree`) for the unit sanity check


## Other Projects

Current list of other [GRI FOSS Projects](https://gitlab.com/geosol-foss/python/gri-convolve/-/blob/main/.docs_other_projects.md) we are building and maintaining.

## License

MIT License. See [LICENSE](https://gitlab.com/geosol-foss/python/gri-convolve/-/blob/main/LICENSE) for details.
