openghg_inversions.basis.algorithms#

Public algorithms for computing basis functions.

This package re-exports weighted and quadtree algorithms alongside constrained basis generation, class composition, allocation, and tuple-safe class masking helpers. It also exposes generic grid-partition contracts, geometry adapters, acceptance policies, and split steps for composing group-local strategies.

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

Class-local strategy derived from recursive weighted bucket splitting.

This applies the existing bucket splitter independently to the masked weights for one class. It recursively splits rectangles along the longer axis until each rectangle is below a searched threshold. It is not a compatibility implementation of the legacy weighted land/sea pipeline, which optimizes the bucket layout before applying the land/sea split. New the default region-constrained strategy is GreedySplitStrategy composed with AxisParallelSplitStep instead.

Variables:: max_iter (int) – Maximum number of threshold-search iterations.

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 GreedySplitStrategy.

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._partition.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#

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

Bases: Protocol

Policy protocol for optional post-construction region consolidation.

class openghg_inversions.basis.algorithms.ConnectedBinaryPartitionStep(split_step: PartitionStep, connectivity: int = 1)#

Bases: object

Repair a provisional binary split into two connected child partitions.

This opt-in wrapper preserves binary arity when disconnected cut fragments can be reassigned safely. It labels the connected components on both sides, retains one primary component on each original side, and moves every secondary component to the opposite side. Candidates are valid only when both resulting children are connected.

Valid candidates are selected deterministically by minimum moved fitting weight, then minimum absolute child-weight imbalance, then row-major component order. If the wrapped step does not return a valid binary partition, or no connected binary reassignment exists, the result falls back to the same multi-child component decomposition as ConnectedComponentPartitionStep.

Variables:

split_step (openghg_inversions.basis.algorithms._partition.PartitionStep) – Partition step whose provisional binary children should be repaired.
connectivity (int) – Two-dimensional neighbourhood definition. 1 uses edge-sharing (four-neighbour) connectivity and 2 additionally includes corner-sharing (eight-neighbour) connectivity.

connectivity: int = 1#

split_step: PartitionStep#

class openghg_inversions.basis.algorithms.ConnectedComponentPartitionStep(split_step: PartitionStep, connectivity: int = 1)#

Bases: object

Make every child from another partition step spatially connected.

This wrapper is useful for partition steps such as InertialSplitStep, whose one-dimensional projection can assign spatially disconnected cells to the same child. Each proposed child is decomposed into deterministic connected components before the greedy orchestrator accepts it.

Variables:

split_step (openghg_inversions.basis.algorithms._partition.PartitionStep) – Partition step whose children should be made connected.
connectivity (int) – Two-dimensional neighbourhood definition. 1 uses edge-sharing (four-neighbour) connectivity and 2 additionally includes corner-sharing (eight-neighbour) connectivity.

connectivity: int = 1#

split_step: PartitionStep#

class openghg_inversions.basis.algorithms.ConnectedComponentSplitStrategy(split_strategy: SplitStrategy, connectivity: int = 1)#

Bases: object

Allocate and split each connected piece of a class independently.

A disconnected class requires at least one label per connected component. If target_regions is below that geographic minimum, the effective target is raised rather than assigning one label to disconnected cells. Targets above the minimum are allocated across components by weight, with the existing area fallback for all-zero weights.

Variables:

split_strategy (openghg_inversions.basis.algorithms._constrained.SplitStrategy) – Class-local strategy applied separately to each connected component.
connectivity (int) – Two-dimensional neighbourhood definition. 1 uses edge-sharing (four-neighbour) connectivity and 2 additionally includes corner-sharing (eight-neighbour) connectivity.
diagnostics (list[dict[str, int]]) – Per-call requested, minimum, effective, and actual region counts.

connectivity: int = 1#

diagnostics: list[dict[str, int]]#

split_strategy: SplitStrategy#

class openghg_inversions.basis.algorithms.ContrastProximityComponentConsolidation(contribution: xr.DataArray | npt.ArrayLike, cell_weight: xr.DataArray | npt.ArrayLike, geometry: LatLonGridGeometry, max_merge_distance_km: float, max_merge_delta_eig: float | None = None, max_merge_lambda: float | None = None, contrast_tau: float | None = None, contrast_sigma_design: float | None = None, contrast_s_diag: xr.DataArray | npt.ArrayLike | None = None, source_classes: xr.DataArray | npt.ArrayLike | None = None, inactive_component_policy: InactiveComponentPolicy = 'keep', connectivity: int = 1, min_regions: int | None = None, spatial_dims: tuple[Hashable, Hashable] | None = None)#

Bases: object

Merge weak, nearby regions created solely by disconnected components.

The input labels must already satisfy the requested strict connectivity. For each source/class group, the largest connected component is treated as the primary component. A secondary component is eligible for consolidation only when it remains one whole basis region. This deliberately excludes one-cell regions created by useful refinement inside the primary component.

Candidate pairs never cross the supplied region_classes or source_classes. Positive-mass candidates are accepted only when every configured reverse-split contrast score is at or below its merge threshold. A zero-mass endpoint is either kept or merged to its nearest eligible neighbour according to inactive_component_policy.

Variables:

contribution (xr.DataArray | npt.ArrayLike) – Fixed design contribution array with one or more design-observation dimensions followed by the two spatial dimensions.
cell_weight (xr.DataArray | npt.ArrayLike) – Non-negative prior mass field used by contrast scoring.
geometry (LatLonGridGeometry) – Latitude/longitude geometry aligned to the label grid.
max_merge_distance_km (float) – Maximum nearest-cell separation for a candidate.
max_merge_delta_eig (float | None) – Optional maximum reverse-split delta_eig.
max_merge_lambda (float | None) – Optional maximum reverse-split lambda.
contrast_tau (float | None) – Prior standard deviation of the split contrast coefficient. None uses the uncalibrated value 1.
contrast_sigma_design (float | None) – Optional scalar design standard deviation.
contrast_s_diag (xr.DataArray | npt.ArrayLike | None) – Optional diagonal design covariance.
source_classes (xr.DataArray | npt.ArrayLike | None) – Optional source field. When omitted, all cells are treated as belonging to one source.
inactive_component_policy (InactiveComponentPolicy) – "keep" or "merge_nearest".
connectivity (int) – 1 for four-neighbour or 2 for eight-neighbour input components.
min_regions (int | None) – Optional global region-count floor.
spatial_dims (tuple[Hashable, Hashable] | None) – Optional contribution spatial dimensions.
diagnostics (list[dict[str, Any]]) – Per-call JSON-compatible consolidation diagnostics.

cell_weight: xr.DataArray | npt.ArrayLike#

connectivity: int = 1#

contrast_s_diag: xr.DataArray | npt.ArrayLike | None = None#

contrast_sigma_design: float | None = None#

contrast_tau: float | None = None#

contribution: xr.DataArray | npt.ArrayLike#

diagnostics: list[dict[str, Any]]#

geometry: LatLonGridGeometry#

inactive_component_policy: InactiveComponentPolicy = 'keep'#

max_merge_delta_eig: float | None = None#

max_merge_distance_km: float#

max_merge_lambda: float | None = None#

min_regions: int | None = None#

source_classes: xr.DataArray | npt.ArrayLike | None = None#

spatial_dims: tuple[Hashable, Hashable] | 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.GreedySplitStrategy(split_step: PartitionStep, split_acceptance: SplitAcceptancePolicy | TargetSplitAcceptancePolicy | None = None)#

Bases: object

Adapt generic greedy partitioning to one class mask.

The strategy converts the selected cells to grid nodes, applies the explicitly supplied partition step through greedy_partitioning(), and converts the result to dense positive class-local labels. It contains no implicit choice of split geometry or algorithm.

Variables:

split_step (openghg_inversions.basis.algorithms._partition.PartitionStep) – Partition step used to propose children for the selected highest-priority partition.
split_acceptance (openghg_inversions.basis.algorithms._partition.SplitAcceptancePolicy | openghg_inversions.basis.algorithms._partition.TargetSplitAcceptancePolicy | None) – Optional policy that may reject a valid proposed split and freeze its parent.

split_acceptance: SplitAcceptancePolicy | TargetSplitAcceptancePolicy | None = None#

split_step: PartitionStep#

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._partition.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, min_child_target_weight_share: float = 0.0)#

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.

By default, every child is subject to the eccentricity limit. Setting min_child_target_weight_share exempts children whose weight is below that share of one class/source-local equal-weight target region, weights.sum() / target_regions. This target-aware exception is used only through accept_split(); the conservative three-argument call remains strict.

The exception affects split acceptance only. It does not reconnect, freeze, prune, or marginalize an exempt child after the split is accepted.

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

Return true when every materially weighted child meets the limit.

weights is the class/source-local field passed to greedy partitioning, so weights.sum() / target_regions is the equal-weight target region weight. Children strictly below min_child_target_weight_share times that reference weight are exempt from the eccentricity veto.

If the total weight is zero, cell counts provide the same direct-call fallback used by MinChildTargetWeightShare; the default greedy strategy already converts all-zero classes to an area surrogate before policies are evaluated. If every child is below the materiality threshold, the strict guard is retained rather than accepting vacuously.

geometry: SplitGeometry | None = None#

max_child_pca_eccentricity: float#

min_child_target_weight_share: float = 0.0#

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 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 – Valid child partitions proposed by a PartitionStep.
weights – Non-negative weight field aligned to the source grid.
target_regions – Requested 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 either input is not two-dimensional, the dimension names differ, weights are invalid, or the requested allocation is impossible.
xarray.AlignmentError – If the inputs do not describe physically compatible spatial grids after transposition.

openghg_inversions.basis.algorithms.combine_inner_outer_region_classes(inner_mask: DataArray, inner_classes: DataArray, outer_classes: DataArray, *, unmapped_values: Iterable[Hashable] = (), name: str = 'region_classes') → DataArray#

Select and tag aligned inner- and outer-domain region classes.

Parameters:

inner_mask – Two-dimensional Boolean field selecting cells from inner_classes. False cells select from outer_classes.
inner_classes – Two-dimensional class field for inner-domain cells.
outer_classes – Two-dimensional class field for outer-domain cells.
unmapped_values – Selected class values that should leave the output cell unmapped. Selected null values are always unmapped.
name – Name for the returned DataArray.

Returns:

Object-valued DataArray with the same dimensions and coordinates as inner_mask. Mapped cells contain ("inner", value) or ("outer", value) tuples. Unmapped cells contain NaN.

Raises:

ValueError – If any input is not two-dimensional, dimension-name sets differ, inner_mask is not Boolean, or a selected class value is not hashable.
xarray.AlignmentError – If, after transposition to inner_mask dimension order, an input does not describe a physically compatible spatial grid.

Notes

Values on the unselected side do not affect the result, including null values. Domain tags prevent equal inner and outer class values from colliding when passed to region_constrained_basis().

This source-neutral composition helper advances issue #449.

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.greedy_partitioning(init_partition: list[list[tuple[int, int]]], target_regions: int, weights: ndarray, *, split_step: PartitionStep, split_acceptance: SplitAcceptancePolicy | TargetSplitAcceptancePolicy | None = None) → list[list[tuple[int, int]]]#

Apply a partition step greedily until a target count is reached.

Parameters:

init_partition – Initial partitions to refine.
target_regions – Refinement ceiling for the number of output partitions. Must be at least one. The function never coarsens an initial partition whose count already exceeds this value.
weights – Two-dimensional non-negative weight field used for partition ranking and passed unchanged to split_step and split_acceptance. Every node must be a valid (row, column) index into this array.
split_step – Callable that proposes child partitions for one selected parent. Fewer than two returned children mark the parent as unsplittable.
split_acceptance – Optional policy applied to a valid multi-child proposal. Rejected splits freeze the selected parent.

Returns:

Output partitions. Empty initial partitions are ignored. The result may contain fewer than target_regions entries when no active partition can be split, an accepted proposal would overshoot the target, or acceptance policies reject the remaining candidates.

Raises:

ValueError – If target_regions is less than one, or a multi-child proposal contains an empty child, duplicate or overlapping nodes, nodes outside its parent, or does not exactly cover its parent.

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. Its region_class_layers attribute records the string form of each layer name in mapping insertion order.

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, after transposition to the first layer’s dimension order, any layer does not describe a physically compatible spatial grid.

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. Use region_class_mask() rather than raw xarray equality when selecting a tuple label.

openghg_inversions.basis.algorithms.normalize_spatial_grid(reference: DataArray, candidate: DataArray, *, reference_name: str = 'reference', candidate_name: str = 'candidate') → DataArray#

Normalize a physically compatible two-dimensional field to a reference grid.

Parameters:

reference – Two-dimensional field whose dimension order, coordinate values, and nonconflicting metadata define the output grid.
candidate – Two-dimensional field to validate and normalize.
reference_name – Name used for reference in validation errors.
candidate_name – Name used for candidate in validation errors.

Returns:

candidate transposed to the reference dimension order and assigned the reference grid coordinates, with nonconflicting coordinate metadata retained from both fields.

Raises:

ValueError – If either field is not two-dimensional or their dimension names differ.
xarray.AlignmentError – If grid coordinates, grid-defining metadata, or CRS definitions are physically incompatible.

Notes

This aligns to an arbitrary reference grid. OpenGHG’s openghg.util.align_lat_lon instead canonicalizes one field against a named OpenGHG domain and does not validate arbitrary curvilinear grids, units, or CRS metadata.

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_class_mask(region_classes: DataArray, class_value: Hashable, *, name: str = 'region_class_mask') → DataArray#

Select one scalar or tuple-valued region class reliably.

Parameters:

region_classes – Region-class values to compare with class_value. Dimensions and all attached coordinates are preserved on the returned mask.
class_value – Scalar or tuple-valued class label to select.
name – Name for the returned Boolean DataArray.

Returns:

Boolean DataArray that is true exactly where region_classes contains class_value.

Notes

Raw expressions such as region_classes == ("land", "inner") are unreliable for object arrays because NumPy may treat the tuple as an array-like operand and broadcast its elements instead of comparing the tuple as one value. This helper performs elementwise tuple comparison through the same path used by constrained allocation and splitting. Object values are materialized eagerly; dask chunks are not preserved in the returned mask.

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, component_consolidation: ComponentConsolidationPolicy | 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 an explicit GreedySplitStrategy using AxisParallelSplitStep.
component_consolidation – Optional policy applied to the globally relabelled basis after class-local construction. Policies that deliberately combine disconnected components must preserve class boundaries and report that strict connectivity no longer holds.
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.

Raises:

ValueError – If either input is not two-dimensional, the dimension names differ, weights are invalid, or the requested allocation is impossible, or a split strategy returns labels with the wrong shape or dtype, non-positive labels inside its class mask, or nonzero labels outside its class mask.
xarray.AlignmentError – If the inputs do not describe physically compatible spatial grids after transposition.

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