openghg_inversions.basis.algorithms#

Algorithms for computing basis functions.

class openghg_inversions.basis.algorithms.AllSplitAcceptancePolicies(*policies: SplitAcceptancePolicy | TargetSplitAcceptancePolicy)#

Bases: object

Accept a split only when every policy accepts it.

accept_split(parent: list[tuple[int, int]], children: list[list[tuple[int, int]]], weights: ndarray, target_regions: int) → bool#: Return true when all component policies accept the split.

policies: tuple[SplitAcceptancePolicy | TargetSplitAcceptancePolicy, ...]#

class openghg_inversions.basis.algorithms.AxisAlignedWeightedSplitStrategy(max_iter: int = 32)#

Bases: object

Compatibility strategy based on the existing recursive weighted basis.

This keeps the current weighted basis shape available for comparison: recursively split rectangles along the longer axis until each rectangle is below a searched threshold. New constrained code defaults to GreedyAxisParallelSplitStrategy instead.

max_iter: int = 32#

class openghg_inversions.basis.algorithms.AxisParallelSplitStep(balanced: bool = True, clean_splits: bool = False, geometry: SplitGeometry | None = None)#

Bases: object

Split one partition along an axis-parallel line.

This is a cleaned-up version of the prototype’s axis-parallel split step. Greedy orchestration is handled separately by GreedyAxisParallelSplitStrategy.

Variables:

balanced (bool) – If true, choose the weighted long axis and split near half total node weight. If false, choose the geometric long axis and split by cell count.
clean_splits (bool) – If true, keep all cells with the same selected-axis coordinate on the same side of the split.
geometry (openghg_inversions.basis.algorithms._constrained.SplitGeometry | None) – Optional geometry used to choose the split axis. The split itself remains a row- or column-aligned cut.

balanced: bool = True#

clean_splits: bool = False#

geometry: SplitGeometry | None = None#

Bases: object

Accept proposed binary splits using an observation-space contrast score.

contribution is a design sensitivity/contribution array with at least one design-observation dimension and native spatial cell dimensions. It is combined with cell_weight to compute h_A = sum_i H_ti * mu_i and h_B = sum_i H_ti * mu_i. If cell_weight is omitted, the class-local weights passed by the greedy splitter are used.

contrast_tau is the prior standard deviation of the new split contrast coefficient, not observation noise. If omitted, tau=1 is used and the score is marked uncalibrated. contrast_sigma_design and contrast_s_diag describe a fixed design covariance in the contribution row space; observed mole-fraction values must not be used here.

If both thresholds are omitted, diagnostics are computable through score_split() but __call__() accepts all valid binary splits.

accepts(score: SplitContrastScore) → bool#: Return true when score satisfies all configured thresholds.

cell_weight: DataArray | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None = None#

contrast_s_diag: DataArray | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None = None#

contrast_sigma_design: float | None = None#

contrast_tau: float | None = None#

contribution: DataArray | _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]#

min_contrast_delta_eig: float | None = None#

min_contrast_lambda: float | None = None#

score_split(child_a: list[tuple[int, int]], child_b: list[tuple[int, int]], *, fallback_cell_weight: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None = None) → SplitContrastScore#: Return contrast-score diagnostics for two child partitions.

spatial_dims: tuple[Hashable, Hashable] | None = None#

class openghg_inversions.basis.algorithms.GreedyAxisParallelSplitStrategy(balanced: bool = True, clean_splits: bool = False, split_step: PartitionStep | None = None, split_acceptance: SplitAcceptancePolicy | TargetSplitAcceptancePolicy | None = None)#

Bases: object

Class-local greedy repeated-bisection strategy.

This is a cleaned-up version of the prototype’s axis-parallel partitioning algorithm. It repeatedly splits the highest-weight current part until the requested class-local region count is reached or no splittable parts remain.

balanced: bool = True#

clean_splits: bool = False#

split_acceptance: SplitAcceptancePolicy | TargetSplitAcceptancePolicy | None = None#

split_step: PartitionStep | None = None#

class openghg_inversions.basis.algorithms.InertialSplitStep(balanced: bool = True, geometry: SplitGeometry | None = None)#

Bases: object

Experimental split step using a weighted principal inertial axis.

The split projects partition cells onto the principal axis of their weighted covariance, then cuts that one-dimensional ordering by weight or by count. This lets diagonal, rotated, or strongly anisotropic high-gradient structures split along their natural orientation instead of being forced through row/column cuts. The greedy class-local orchestrator still invokes this step independently inside each region class, so labels keep the same region-constrained boundary guarantees as axis-parallel splitting.

By default the covariance uses grid-index coordinates. Pass geometry=LatLonGridGeometry.from_dataarray(...) to use local physical north-south and east-west metre offsets for each selected partition. Degenerate covariance, tied projections at the selected cut, and other numerically unstable cases fall back to an axis-parallel split.

Variables:

balanced (bool) – If true, split the inertial projection near half total node weight. If false, split by cell count.
geometry (openghg_inversions.basis.algorithms._constrained.SplitGeometry | None) – Optional geometry used for the covariance and projection. The fallback split uses the same geometry for axis choice.

balanced: bool = True#

geometry: SplitGeometry | None = None#

class openghg_inversions.basis.algorithms.LatLonGridGeometry(latitudes: ndarray[tuple[Any, ...], dtype[float64]], longitudes: ndarray[tuple[Any, ...], dtype[float64]], earth_radius_m: float = 6371008.8)#

Bases: object

Local tangent-plane geometry for latitude/longitude grids.

Coordinates are computed per partition in metres using a local equirectangular approximation centered on the weighted latitude/longitude of the selected nodes. The returned coordinate columns are local north-south and east-west metre offsets, matching grid axes 0 and 1 for row/column split decisions.

Variables:

latitudes (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Finite two-dimensional latitude coordinate grid in degrees, aligned to grid node (row, col) indexing.
longitudes (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Finite two-dimensional longitude coordinate grid in degrees, with the same shape and node alignment as latitudes.
earth_radius_m (float) – Earth radius used for converting angular differences to metres. Must be positive and finite.

coordinates(nodes: list[tuple[int, int]], node_weights: ndarray[tuple[Any, ...], dtype[float64]] | None = None) → ndarray[tuple[Any, ...], dtype[float64]] | None#

Return local tangent-plane coordinates for nodes in metres.

Parameters:

nodes – Grid nodes in the selected partition.
node_weights – Optional non-negative weights for the same nodes. These weights set the local projection center. Equal weights are used when weights are omitted or all zero.

Returns:

A finite (nnode, 2) array of local north-south and east-west metre offsets. The local center is the weighted mean latitude and circular weighted mean longitude for nodes, so partitions near the antimeridian use the shorter wrapped longitude difference. Empty nodes returns an empty coordinate array. Invalid coordinates, out-of-bounds nodes, or invalid node_weights return None so callers can fall back to row/column index coordinates.

earth_radius_m: float = 6371008.8#

classmethod from_dataarray(data: DataArray, *, lat_name: str = 'lat', lon_name: str = 'lon', earth_radius_m: float = 6371008.8) → LatLonGridGeometry#

Create geometry from latitude and longitude coordinates.

Parameters:

data – Two-dimensional grid with latitude and longitude coordinates, ordered as (lat_name, lon_name) so node axes align with the returned north-south and east-west metre offsets.
lat_name – Name of the latitude coordinate and first dimension.
lon_name – Name of the longitude coordinate and second dimension.
earth_radius_m – Earth radius used for metre scaling.

Returns:

Geometry aligned to data.

Raises:

ValueError – If dimensions are not ordered as (lat_name, lon_name) or coordinates cannot be broadcast to the data grid.

latitudes: ndarray[tuple[Any, ...], dtype[float64]]#

longitudes: ndarray[tuple[Any, ...], dtype[float64]]#

class openghg_inversions.basis.algorithms.MaxChildPCAEccentricity(max_child_pca_eccentricity: float, geometry: SplitGeometry | None = None, tolerance: float = 1e-12)#

Bases: object

Reject splits that create child partitions above a PCA eccentricity limit.

The eccentricity is computed from each child partition’s unweighted node coordinates. If geometry is supplied, its physical coordinates are used; otherwise row/column index coordinates are used. Single-cell children have eccentricity 1 because they have no resolvable long axis. Multi-cell rank-one children have infinite eccentricity and are rejected by any finite threshold.

geometry: SplitGeometry | None = None#

max_child_pca_eccentricity: float#

tolerance: float = 1e-12#

class openghg_inversions.basis.algorithms.MinChildTargetWeightShare(min_child_target_weight_share: float)#

Bases: object

Reject splits whose lightest child is below an equal-target share.

min_child_target_weight_share is compared with min(child_weight) / (weights.sum() / target_regions) for the class/source-local weights being partitioned. This policy stops creation of low-weight basis regions relative to the requested equal-weight target; it is not a parent-relative split-balance guard.

accept_split(parent: list[tuple[int, int]], children: list[list[tuple[int, int]]], weights: ndarray, target_regions: int) → bool#

Return true when every child is large enough to become a region.

weights is the class/source-local field passed to greedy partitioning, so weights.sum() / target_regions is the equal-weight target region weight. If the total weight is zero, fall back to cell-count shares for direct policy use; the default greedy strategy already converts all-zero classes to an area surrogate before policies are evaluated.

min_child_target_weight_share: float#

class openghg_inversions.basis.algorithms.MinChildWeightShare(min_child_weight_share: float)#

Bases: object

Reject splits whose lightest child is below a parent-weight share.

This is a split-balance guard. It compares children with their current parent partition, not with the total class/source weight being partitioned.

min_child_weight_share: float#

class openghg_inversions.basis.algorithms.PartitionStep(*args, **kwargs)#

Bases: Protocol

Strategy protocol for splitting one partition into child partitions.

class openghg_inversions.basis.algorithms.SplitAcceptancePolicy(*args, **kwargs)#

Bases: Protocol

Policy protocol for accepting proposed child partitions.

class openghg_inversions.basis.algorithms.SplitContrastScore(contrast: ndarray[tuple[Any, ...], dtype[float64]], lambda_value: float, delta_dfs: float, delta_eig: float, mu_a: float, mu_b: float, tau: float, uncalibrated: bool)#

Bases: object

Diagnostics for one mass-preserving split contrast.

Variables:

contrast (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – The contrast column f_ab in design-observation space.
lambda_value (float) – tau**2 * f_ab.T @ S^{-1} @ f_ab.
delta_dfs (float) – Incremental DFS proxy lambda / (1 + lambda).
delta_eig (float) – Incremental EIG proxy 0.5 * log(1 + lambda).
mu_a (float) – Prior mass in child A.
mu_b (float) – Prior mass in child B.
tau (float) – Prior standard deviation of the split contrast coefficient delta = alpha_A - alpha_B.
uncalibrated (bool) – True when the score used the default tau=1 or identity covariance. Such scores are useful for ranking/debugging but not calibrated expected information gain.

contrast: ndarray[tuple[Any, ...], dtype[float64]]#

delta_dfs: float#

delta_eig: float#

lambda_value: float#

mu_a: float#

mu_b: float#

tau: float#

uncalibrated: bool#

class openghg_inversions.basis.algorithms.SplitGeometry(*args, **kwargs)#

Bases: Protocol

Geometry protocol for mapping grid nodes into physical coordinates.

coordinates(nodes: list[tuple[int, int]], node_weights: ndarray[tuple[Any, ...], dtype[float64]] | None = None) → ndarray[tuple[Any, ...], dtype[float64]] | None#

Return physical coordinates for grid nodes.

Parameters:

nodes – Grid nodes in the partition being split.
node_weights – Optional non-negative weights for the same nodes.

Returns:

A finite (nnode, 2) coordinate array whose first column is aligned with grid axis 0 and second column is aligned with grid axis 1. Return None when physical coordinates are unavailable and index-space fallback should be used.

class openghg_inversions.basis.algorithms.SplitStrategy(*args, **kwargs)#

Bases: Protocol

Strategy protocol for class-local basis splitting.

class openghg_inversions.basis.algorithms.TargetSplitAcceptancePolicy(*args, **kwargs)#

Bases: SplitAcceptancePolicy, Protocol

Policy protocol for accepting splits using the class-local target count.

accept_split(parent: list[tuple[int, int]], children: list[list[tuple[int, int]]], weights: ndarray, target_regions: int) → bool#

Return true when proposed children should be accepted.

Parameters:

parent – Parent partition selected by greedy orchestration.
children – Non-empty child partitions proposed by a PartitionStep.
weights – Non-negative class/source-local weight field aligned to the source grid.
target_regions – Requested class/source-local upper target count.

Returns:

True if greedy orchestration should replace parent with children. False freezes parent as a completed partition.

openghg_inversions.basis.algorithms.allocate_nbasis_by_class(weights: DataArray, region_classes: DataArray, nbasis: int | Mapping[Hashable, int], *, allocation: Literal['weight', 'area'] = 'weight', min_regions_per_class: int = 1, unmapped_values: Iterable[Hashable] = ()) → dict[Hashable, int]#

Allocate class-local region targets for constrained basis generation.

Parameters:

weights – Two-dimensional non-negative weight field.
region_classes – Two-dimensional class field aligned to weights.
nbasis – Total number of regions to distribute, or an explicit mapping from class value to class-local region target.
allocation – Automatic allocation mode. "weight" uses class total weight, falling back to area if all mapped weights are zero. "area" uses mapped cell count.
min_regions_per_class – Minimum automatic allocation for each non-empty mapped class.
unmapped_values – Additional class values to leave unallocated.

Returns:

Mapping from mapped class value to target number of local regions.

Raises:

ValueError – If inputs cannot be aligned, weights are invalid, or the requested allocation is impossible.

openghg_inversions.basis.algorithms.contrast_tau_from_multiplier_cv(multiplier_cv: float, *, approximation: Literal['additive', 'log'] = 'additive') → float#

Return an approximate split-contrast tau from a multiplier CV.

tau is the prior standard deviation of delta = alpha_A - alpha_B. The additive approximation uses sqrt(2) * multiplier_cv. The log approximation uses sqrt(2) * sqrt(log1p(multiplier_cv**2)).

openghg_inversions.basis.algorithms.intersect_region_class_layers(layers: Mapping[Hashable, DataArray], *, unmapped_values: Iterable[Hashable] = (), name: str = 'region_classes') → DataArray#

Intersect aligned region-class layers into composite class labels.

Parameters:

layers – Ordered mapping from layer name to two-dimensional class field. Mapping insertion order defines the order of values in each output class tuple.
unmapped_values – Layer values that should leave the output cell unmapped. Null values in any layer are always unmapped.
name – Name for the returned DataArray.

Returns:

Object-valued DataArray with the same dimensions and coordinates as the first layer. Mapped cells contain tuples of layer values, while cells that are null or explicitly unmapped in any layer contain NaN.

Raises:

ValueError – If no layers are supplied, any layer is not two-dimensional, layer dimension names differ, or a mapped layer value is not hashable.
xarray.AlignmentError – If layer coordinates do not align exactly.

Notes

The tuple labels can be passed directly to region_constrained_basis(). This is the small lattice-style construction needed for layered masks such as land/sea crossed with an inner/outer rectangle.

openghg_inversions.basis.algorithms.quadtree_algorithm(fps: ndarray, nbasis: int, seed: int | None = None) → ndarray#

Given an array and a specified number of basis functions, return basis regions specified by the quadtree algorithm.

Parameters:

fps – array (mean flux times mean footprints) to use to calculate basis regions
nbasis – target number of basis regions
seed – optional random seed to use (for testing or reproducing results)

Returns:

2D numpy array with positive integer values representing basis regions.

openghg_inversions.basis.algorithms.region_constrained_basis(weights: DataArray, region_classes: DataArray, nbasis: int | Mapping[Hashable, int], *, allocation: Literal['weight', 'area'] = 'weight', min_regions_per_class: int = 1, split_strategy: SplitStrategy | None = None, unmapped_values: Iterable[Hashable] = ()) → DataArray#

Generate basis labels independently inside each mask/region class.

Parameters:

weights – Two-dimensional non-negative weight field.
region_classes – Two-dimensional class field on the same grid as weights. Each non-null value is treated as a mapped class unless listed in unmapped_values.
nbasis – Either a total number of basis regions to allocate across classes, or an explicit mapping from class value to class-local region target.
allocation – Automatic allocation mode used when nbasis is an integer. "weight" allocates proportional to class total weight, falling back to area if all class weights are zero. "area" allocates proportional to mapped cell count.
min_regions_per_class – Minimum automatic allocation for each non-empty mapped class. If the requested total is smaller than this minimum requires, a ValueError is raised.
split_strategy – Class-local splitting strategy. Defaults to GreedyAxisParallelSplitStrategy.
unmapped_values – Additional class values to leave as output label 0.

Returns:

xarray.DataArray with the same dimensions and coordinates as weights. Mapped cells receive globally unique positive integer labels; unmapped cells receive 0.

Notes

Labels are guaranteed not to cross class boundaries because each class is split independently and relabelled with a global offset. The default strategy can assign one label to disconnected pieces of the same class if the class mask itself is disconnected; contiguity is not guaranteed by this helper.

Compute the mass-preserving split contrast score.

contribution must have at least one design-observation dimension plus the native spatial cell dimensions. cell_weight supplies the positive prior flux/mass weights mu_i. The score is based only on design sensitivities/contributions and prior mass; observed mole fractions or residuals are not inputs.

openghg_inversions.basis.algorithms.weighted_algorithm(grid: ndarray, bucket: float = 1, nregion: int = 100, tol: int = 1, domain: str = 'EUROPE', country_directory: str | None = None) → ndarray#

Obtain basis function with nregions (for land-sea split).

Parameters:

grid – 2D grid of footprints * flux, or whatever grid you want to split. Could be: population data, spatial distribution of bakeries, you choose!
bucket – Initial bucket value for each basis function region. Defaults to 1
nregion – Number of desired basis function regions Defaults to 100
tol – Tolerance to find number of basis function regions. i.e. optimizes nregions to +/- tol Defaults to 1
domain – Domain across which to calculate basis functions.
country_directory – Directory containing land-sea files. If None, will use default files.

Returns:

2D basis function array

Return type:

basis_function