Bins, tapers, fields and geometry

AbstractBinEdges{T} <: AbstractVector{T}

Supertype for all custom, high-performance bin edge collections in StructureFunctions.jl.

Why AbstractBinEdges Exists

In structure function calculations over large datasets, the spatial separation distance $r$ for each of the $O(N^2)$ point pairs must be mapped to its corresponding distance bin (index). Using a standard sorted vector of bin edges requires a binary search (searchsortedfirst), which has $O(\log B)$ complexity where $B$ is the number of bins.

For large $N$, this binary search becomes the dominant CPU bottleneck, causing high branch mispredictions and cache misses. Subtypes of AbstractBinEdges bypass the standard binary search by implementing custom Base.searchsortedfirst overrides that execute in $O(1)$ time:

• LinearBinEdges utilizes Fused Multiply-Add (FMA) arithmetic for uniformly-spaced bins.
• LogBinEdges maps physical queries via log(q) then the same FMA path on the log-space grid.

Wrapping standard arrays in these subtypes makes digitize a constant-time lookup.

Bin edges and backends

• Plain AbstractVector bin edges are valid everywhere but use generic searchsortedfirst / binary search.
• Pass LinearBinEdges, LogBinEdges, or InfPaddedBinEdges for O(1) CPU digitize and matching GPU tiled kernels (see StructureFunctionsKernelAbstractionsExt).
• InfPaddedBinEdges adds implicit catch-all under/overflow bins; do not vcat(-Inf, …, Inf) manually.
• Classic serial/threaded pair-loop paths normalize once via BinEdges; single-pass APIs expect callers to choose edge types explicitly (no auto-wrap on GPU).

AbstractTaper

A weight applied before a transform: on a lag-space autocovariance, as a function of the lag's length; on a spherical harmonic series, as a function of the degree. NoTaper() leaves every term as it is, Bartlett() falls linearly to zero at the largest lag or degree, GaussianTaper(σ) weights a lag of length r by exp(-r²/2σ²) and a degree l by exp(-l(l+1)σ²/2), with σ in the lag's or the sphere's own unit. A taper trades resolution for variance, or a hard bin for a positive kernel.

The taper falling linearly to zero at the largest lag or degree; see AbstractTaper.

BinEdges(edges::AbstractVector{T})

Generic wrapper for arbitrary sorted vectors of bin edges. Bypasses range-specific optimizations but conforms to the AbstractBinEdges interface.

Behavior

• If constructed with an AbstractRange (e.g. StepRange or StepRangeLen), it automatically promotes and returns a LinearBinEdges wrapper to enable O(1) FMA indexing.
• Otherwise, it wraps the vector and delegates to standard $O(\log N)$ binary search methods.

BinEdges(edges)

Normalize flat edge input to AbstractBinEdges for hot-loop digitize:

• existing AbstractBinEdges → unchanged
• AbstractRange → LinearBinEdges
• other AbstractVector → wrapped via the default BinEdges struct constructor (generic binary search)

GaussianTaper(σ)

The Gaussian taper exp(-r²/2σ²) on a lag of length r, exp(-l(l+1)σ²/2) on a degree l, or exp(-k²σ²/2) on a mode of wavenumber k, with σ a length; see AbstractTaper.

HarmonicNodes(separations, lmax; taper = NoTaper())
HarmonicNodes(n::Integer, lmax; taper = NoTaper())

The "bins" of a kernel-binned pair statistic on a sphere: the central angles separations (radians) at which it is reported, the degree lmax its harmonic series is truncated at, and the taper on that series. Together they are the kernel, K(γ, β) = (1/16π²) Σ_{l ≤ lmax} (2l+1) b_l d^l(cos γ) d^l(cos β), which replaces a hard bin around β; it narrows as π/lmax and tends to a delta in cos γ.

weights are quadrature weights in μ = cos β, ∫_{-1}^{1} f dμ ≈ Σ_k w_k f(μ_k): Gauss–Legendre for the n-node form, whose nodes are the Gauss–Legendre points of μ, and the midpoint rule on the nodes' own cells for given separations. They are what inverts the statistic back to a spectrum.

Not a vector of edges: one value is reported per node, so a result built on these has as many values as nodes.

InfPaddedBinEdges(edges::AbstractVector{T})

Wrapper that implicitly prepends $-\infty$ (or typemin(T)) and appends $+\infty$ (or typemax(T)) to an existing bin edge collection.

Why InfPaddedBinEdges Exists

Structure function distance bins are defined as half-open intervals $(r_i, r_{i+1}]$. When mapping a distance $r$ to a bin, any query value $r < \text{first}(edges)$ or $r > \text{last}(edges)$ is out-of-bounds. InfPaddedBinEdges embeds the infinite endpoints, so an inner loop needs no out-of-bounds branch:

• The first element is treated as typemin(T) ($-\infty$).
• The last element is treated as typemax(T) ($+\infty$).

This guarantees that every valid positive separation distance maps to a valid index without allocating actual padding elements in memory or copying the array.

Prevention of Double Padding

The constructor checks if the input array already has infinite endpoints. If they exist, it trims them before wrapping to prevent nested padding (e.g. $[-\infty, -\infty, ...]$).

LinearBinEdges(edges::AbstractRange{T})

High-performance wrapper for uniformly-spaced ranges (linear spacing).

Mathematical Theory

A standard binary search takes $O(\log B)$ steps. With a uniformly spaced range of bin edges $v_i = v_1 + (i-1)\delta$, searchsortedfirst(v, x) is the smallest index $i$ with $v_i \ge x$:

\[i^*(x) = \min\{ i : v_1 + (i-1)\delta \ge x \} = \left\lfloor \frac{x - v_1}{\delta} \right\rfloor + 1 = \lceil \frac{x - v_1}{\delta} + 1 \rceil .\]

The discrete operator is ceiling / floor+1, not round (which answers a different question).

Precompute inv_step = 1/δ and at query time one FMA gives $t = (x-v_1)/\delta$: t = muladd(x, inv_step, -v_1 * inv_step) then idx = floor(Int, t) + 1.

Float boundary correction

t and reconstructed edges are inexact. After clamping idx, compare the reconstructed left edge u_idx = v_1 + (idx-1)δ; if u_idx < x, return idx + 1. One comparison fixes the at-most-one-bin FP error without a second downward correction.

Performance

Bypasses the twice-precision arithmetic of Julia's StepRangeLen search, so a lookup is one FMA and one comparison.

LogBinEdges(edges::AbstractVector{T})
LogBinEdges_from_log_edges(log_edges)

Log-spaced (geometric) bin edges: uniform grid in log-space.

Digitize semantics

log_edges is the authoritative grid. A physical query q > 0 maps to searchsortedfirst(LinearBinEdges(log_edges), log(q)). See the documentation's "Binning Internals" page.

LogBinEdges from a physical vector builds the same log grid from range(log(first), log(last); length). LogBinEdges_from_log_edges accepts the log grid directly.

Performance

One log(q) plus the constant-time FMA digitize on the log grid.

ModeBinEdges(edges, schedule)

Distance bins applied to the lags of a mode grid. The edges are ordinary bin edges, but every pair reaches them through the periodic kernel of the mode set schedule describes, so a bin's sum and count are kernel-weighted and the histogram is soft-binned. Carried by the results of the non-uniform FFT route, so they cannot be mistaken for pair counts.

The taper that leaves every term as it is; see AbstractTaper.

LogBinEdges_from_log_edges(log_edges)

Build log-uniform bins from log-space edges (u with physical edge exp(u)).

gauss_legendre(n) -> (nodes, weights)

The n Gauss–Legendre nodes on (-1, 1), ascending, and their weights, by Newton's method on P_n from the Tricomi estimate.

The outer bins of InfPaddedBinEdges are unbounded, so they have no representative separation. Take midpoints of the bounded interior, edges.edges.

midpoints(edges::LinearBinEdges) -> LinearBinEdges
midpoints(edges::LogBinEdges) -> LogBinEdges

Midpoints of LinearBinEdges / LogBinEdges.

n_histogram_bins(edges::AbstractVector) -> Int

Number of histogram bins for flat edges (length == N + 1 → N bins).

Bins covered by the plan (digitize results in 1:n_bins are in range).

physical_edges_vector(bins::LogBinEdges) -> Vector

Materialize physical bin edges exp(log_edges[i]) for display or generic consumers. Not used on the digitize hot path.

Several fields at one set of points, differenced together.

A structure function is a statement about one pair and any number of quantities carried at its two ends. Velocity is a vector field, parallel-transported by the geometry before differencing; a tracer, a vorticity, a temperature is a scalar field, differenced as it stands. Both are the same pair sweep, so they travel together in one Fields.

FieldIncrement(vectors, scalars)

One pair's increment across every field: each vector field already transported into the pair's common frame, each scalar field already differenced.

An operator reads the fields it names. A field of one vector field and no scalars does not produce one of these — its increment is the plain SVector every existing operator already takes, so the single-field path is unchanged down to the instruction.

Fields(; vectors = (), scalars = ())

Several quantities sampled at one set of points: vectors are parallel-transported before being differenced, scalars are differenced as they stand.

Fields(vectors = (u,))                        # a velocity — what a bare `u` already means
Fields(vectors = (u,), scalars = (θ,))        # velocity and a tracer, for the mixed laws
Fields(vectors = (u, 𝓐u))                     # velocity and its advection
Fields(scalars = (ω, 𝓐ω))                     # a scalar and its advection

Each vector is (D, N) and each scalar is (N,) or (1, N); on a grid, each vector is (D, cells...) and each scalar (cells...), the cells flattened in the order the array stores them. They are packed once, at construction, into the single (V·D + K, N) array the kernels already load — so a pair costs one contiguous read whatever it carries, and the field counts are type parameters, so the kernel specialises on them.

Which is transported and which is not follows from the name and from nothing else: a vector field carries a direction the geometry must rotate, a scalar field does not.

Components in each vector field.

Tracer-like fields, each differenced without transport.

Velocity-like fields, each parallel-transported.

Generic helper functions for the other routines in this package

AbstractTransverseBasisConvention

A rule giving the transverse basis of a pair from its unit separation: transverse_basis(rule, r̂) returns a tuple of unit vectors perpendicular to r̂ and to each other, whose first vector is odd under r̂ ↦ −r̂. Operators odd in the transverse component read that first vector, and its oddness is what makes them read the same from either end of a pair. A new rule is a subtype with that one method.

CanonicalTransverseBasis()

The transverse direction n̂: in 2-D the quarter turn of r̂; in 3-D the turn about ẑ, continued about x̂ where r̂ ∥ ẑ, followed by r̂ × n̂. Every projected operator carries this convention unless constructed with another.

FlatGeometry{D}()

Flat D-dimensional space with the Euclidean ruler: the separation is the straight chord and no transport is needed.

ReferenceAxisTransverseBasis(axis; parallel_tol = 1e-12)

The rule of n̂ about a chosen axis a in 3-D: e₁ = normalize(a × r̂) and e₂ = r̂ × e₁. With a = ẑ it is CanonicalTransverseBasis without the continuation along ẑ: where ‖a × r̂‖² ≤ parallel_tol² it throws ArgumentError. Refused in 2-D, where the transverse line is n̂ and an axis adds nothing.

SphericalDistance(radius)

Great-circle distance on a sphere of the given radius between points given as (longitude, latitude) in radians: radius times Distances.SphericalAngle. The one metric a spherical grid of known radius hands to every route, so separations come out in the radius's unit on each of them.

SphericalGeometry{D}(metric, radius)

Sphere of the given radius. D is the velocity dimension: 2 for horizontal (east, north) velocities, 3 for a thin shell carrying an additional radial component.

D cannot be inferred from the coordinates — a point on a shell has two of them either way — and it changes the answer, because transverse_component_norm2 divides the transverse energy by D - 1.

metric is retained because the coordinates mean whatever it says they mean: it is what fixes their angle unit, which pair_frame needs in order to take a sincos. See unit_position.

coordinate_width(geometry) -> Val{W}

How many numbers locate one point in the form the kernels consume, as a Val so callers can build a statically-sized load. It is the coordinate count, never the velocity dimension: a point on a shell takes three ambient components whether or not the velocity carries a radial one.

See input_coordinate_width for the width a caller supplies and prepare_pair_inputs for the conversion between the two.

digitize(x, bins) -> Int

The index of the bin (bins[i], bins[i + 1]] holding x: 0 below the first edge and length(bins) above the last. A vector x gives one index per element. The bins are right-inclusive.

field_width(geometry) -> Val{W}

How many components of a field the kernels index per point.

Defaults to coordinate_width: in the kernel representation a position and a field are ambient vectors of the same dimension — D apiece in flat space, three apiece on a sphere, where prepare_pair_inputs has expressed both in ambient coordinates.

The geometry's D stays the velocity dimension and is what selects the transverse normalization and whether a radial component is carried, so it cannot be read off a prepared array's leading axis.

flatten_data(vecs::Tuple)

Flatten multi-dimensional arrays in the tuple to vectors. Useful for converting grid data to point lists.

geodesic_frame(p̂, q̂) -> (σ, t̂_A, t̂_B, m̂, ok)

Great-circle frame for the pair (p̂, q̂) of unit position vectors: central angle σ, the longitudinal unit tangent at each endpoint (t̂_A at p̂ toward q̂, t̂_B at q̂ continuing along the geodesic), and the shared transverse unit vector m̂.

m̂ is the same 3-vector at both endpoints: the great-circle normal is perpendicular to every point of that circle, hence tangent to the sphere at each of them, and it is parallel along the geodesic. Both tangents also share one normalizer, since ‖q̂ − (p̂·q̂)p̂‖² = ‖(p̂·q̂)q̂ − p̂‖² = 1 − (p̂·q̂)². So the whole frame costs one sqrt.

σ uses the tangent-half-angle form 2·atan(‖p̂−q̂‖, ‖p̂+q̂‖), accurate for every σ including antipodal. acos(p̂·q̂) loses half the mantissa near σ=0, which is fatal in Float32.

ok is false for coincident (σ=0) and antipodal (σ=π) pairs, where the direction is genuinely undefined — antipodal points are joined by infinitely many great circles, so parallel transport between them is not unique. Callers must skip such pairs: 1/s is otherwise Inf and a single NaN lane poisons an entire atomic accumulator.

geodesic_increments(t̂_A, t̂_B, m̂, u_A, u_B) -> (δu_L, δu_T)

Longitudinal and transverse components of the parallel-transported velocity difference, given the frame from geodesic_frame and the two ambient tangent velocities.

Projecting each endpoint's velocity onto its own geodesic tangent is parallel transport, not an approximation: a geodesic parallel-transports its own tangent, and transport on the sphere is an orientation-preserving isometry, so the transported basis from A arrives at B rotated by exactly the difference of forward azimuths. The transverse term needs a single dot product because m̂ is shared between the endpoints.

input_coordinate_width(geometry) -> Val{W}

How many numbers a caller supplies per point: a longitude and a latitude on a sphere, which prepare_pair_inputs turns into the three-component unit position the kernels index.

Defaults to coordinate_width, so a geometry whose input needs no conversion implements nothing.

local_east_north(lon, lat) -> (Ê, N̂)

Local east and north unit vectors at longitude/latitude in degrees. Together with unit_position as the local up they form a right-handed orthonormal triad (Ê × N̂ = p̂). Only used at ingest: p̂ itself is well defined at the poles, this basis is not.

Return the signed longitudinal magnitude of δu along r_hat. The caller must ensure r_hat is a unit vector.

Return the transverse magnitude of δu, signed relative to the normal vector n̂(r_hat). The caller must ensure r_hat is a unit vector.

Return the transverse (perpendicular) unit vector from x1 to x2, defined as the cross product of the local vertical with the longitudinal unit vector. This naming is opposite to the Lindborg and Cho convention.

n̂(r_hat)

Oriented transverse unit vector for the longitudinal unit vector r_hat.

2D: n̂ = ẑ × r̂ = (−r̂₂, r̂₁), the counterclockwise quarter turn, so (r̂, n̂, ẑ) is right-handed. 3D: n̂ = normalize(ẑ × r̂), the same rule with ẑ = (0,0,1) as the reference axis; along ẑ, where that product vanishes (‖ẑ × r̂‖² ≤ eps), the rule continues with x̂ as the axis, n̂ = normalize(x̂ × r̂), which is −ŷ at r̂ = ẑ and +ŷ at r̂ = −ẑ.

Only operators odd in the transverse component see this sign — ProjectedStructureFunctionType{2,1} and {0,3}. Everything else consumes δu_T² (see transverse_norm2) and is sign-blind. As a convention carried by an operator this rule is CanonicalTransverseBasis; one with a chosen axis is ReferenceAxisTransverseBasis.

pair_delta(geometry, frame, x1, x2, u1, u2) -> δu

Velocity difference expressed in the pair's common frame, given the frame from pair_frame. The per-velocity-field half of pair_increments.

Increments on the sphere from the ambient velocities prepare_pair_inputs produces. Each endpoint's velocity is projected onto its own geodesic tangent and only then differenced — that IS the parallel transport, because transport along a geodesic preserves the angle to it. m̂ is shared between the endpoints, so the transverse term is a single dot product.

Thin shell. The radial component is recovered as u·p̂ at each endpoint and differenced as a scalar: t̂ ⟂ p̂, q̂ and m̂ ⟂ p̂, q̂, so it cannot leak into the tangential increments and needs no transport.

pair_direction(geometry, frame, r) -> r̂

Longitudinal unit vector for one pair. Depends only on the geometry of the pair, never on the velocities, so a caller sweeping many velocity fields over one (i, j) hoists this out of that loop.

On the sphere the frame IS the basis, so the longitudinal direction is ê₁ by construction.

pair_frame(geometry, x1, x2) -> (ok, r, frame)

Separation and the geometry-specific frame data for one pair, touching no velocities.

Split from pair_increments because callers depend on that order: the scalar kernels reject out-of-range pairs before loading u_j at all, and the batch kernels build one frame per (i, j) and reuse it across a whole strip of velocity fields. A single call taking the velocities would defeat both.

Spherical frame from the ambient unit positions

pair_geometry(geometry, x1, x2, u1, u2) -> (ok, r, δu, r̂)

Resolve one pair into a separation r, a velocity difference δu, and a longitudinal unit vector r̂, all in one common frame, so that every structure-function operator — which consumes only δu·r̂ and δu·n̂(r̂) — works unchanged for any geometry. ok == false marks a pair with no defined separation direction; callers must skip it.

For SphericalGeometry, x1/x2 are (lon, lat) in the metric's angle unit and u1/u2 are (east, north[, up]). The result is expressed in the local geodesic frame, so r̂ = ê₁ and δu = (δu_L, δu_T) (plus the radial difference for D = 3). The radial component needs no transport and cannot leak into δu_L/δu_T, because t̂ ⟂ p̂ and m̂ ⟂ p̂, q̂.

This is the one-shot convenience form; hot loops use pair_frame and pair_increments separately so the frame can gate the bin test and be reused across fields.

pair_geometry_for(metric, ::Val{D}) -> geometry

Geometry implied by metric for user-facing dimension D. A distance function alone does not define a direction or a transport rule, so there is deliberately no generic method and an unrecognized metric raises. Add a method here (and a pair_geometry method for the geometry it returns) to support another manifold.

pair_increments(geometry, frame, r, x1, x2, u1, u2) -> (δu, r̂)

Velocity difference and longitudinal unit vector in one common frame, given the frame and separation r from pair_frame. See pair_geometry for what the result means per geometry. Called only for pairs that survive the bin test, so this is where the direction is formed.

pair_invariants(geometry, frame, r, u1, u2) -> (δu_L, ‖δu‖²)

The only two scalars the six isotropic invariants consume: every one of them is built from δu_L, δu_L² and δu_T² = ‖δu‖² − δu_L². A kernel that needs no separation direction calls this and never forms r̂; pair_increments is the form that does.

pair_orientation(geometry, frame) -> Int

Which end of a pair is read first: +1 when the pair as given runs from the lower to the upper end along the first coordinate that separates its points, -1 when it runs the other way, 0 when neither end comes first. On a sphere a pair is read from south to north, and along one parallel from west to east by the shorter way round. Only operators odd in a scalar increment see this sign; every vector quantity is the same read either way.

prepare_coordinates(geometry, x) -> x_kernel

The coordinate half of prepare_pair_inputs, for a caller that has no field to convert alongside — a field of scalar fields only has nothing to transport, but its points still have to reach the form the kernels index.

prepare_pair_inputs(geometry, x, u) -> (x_kernel, u_kernel)

Convert caller coordinates and fields into the form the kernels index, once per point.

Flat geometry passes both through. Spherical geometry maps (lon, lat) to the ambient unit position p̂, and (east, north[, up]) to the ambient velocity u₃ = u_E Ê + u_N N̂ [+ u_r p̂], both (3, N).

Trailing auxiliary axes are preserved. x may carry one slice against many field slices (the fixed-position batch), in which case its single slice supplies the basis for all of them.

remove_nans(x_mat::AbstractMatrix, u_mat::AbstractMatrix)

Remove points (columns) where either the position x_mat or the velocity u_mat contains a NaN. Returns (x_mat_clean, u_mat_clean).

Return the longitudinal (parallel) unit vector from x1 to x2.

transverse_basis(rule::AbstractTransverseBasisConvention, r̂) -> NTuple of unit vectors

The transverse basis of a pair with unit separation r̂ under rule: one vector in two dimensions, two in three, each perpendicular to r̂ and to the others, the first odd under r̂ ↦ −r̂.

transverse_basis_vector(r̂, rule, index = 1)

The index-th vector of transverse_basis(rule, r̂); the first is the signed transverse direction an odd transverse operator reads.

transverse_component(δu, r̂, basis, basis_index=1)

Component of δu along an explicit transverse basis vector. This is the quantity to use for basis-dependent odd/component diagnostics such as T3SF or L2T1SF.

transverse_component_norm2(δu, r̂)

Per-component transverse energy under isotropic transverse averaging, $||δu_t||² / (D - 1)$.

transverse_norm2(δu, r̂)

Squared norm of the full transverse vector $δu - (δu⋅r̂)r̂$. This is the invariant transverse energy used by T2SF and L1T2SF.

unit_position(lon, lat) -> SVector{3}

Unit position vector on the sphere from longitude/latitude in degrees (the Distances.Haversine convention). x toward (0°N, 0°E), z toward the north pole.

unit_position(metric, lon, lat) -> SVector{3}
local_east_north(metric, lon, lat) -> (Ê, N̂)

Ingest helpers that take the angle unit from the metric's own documented convention: Distances.Haversine is degrees, Distances.SphericalAngle and SphericalDistance are radians. Confusing the two silently rescales every separation by a factor of ~57, so the convention is pinned next to the metric that defines it.

Return the vector from x1 to x2.

Return the longitudinal component of δu along r_hat, as a vector. The caller must ensure r_hat is a unit vector.

Return the transverse component of δu (perpendicular to r_hat), as a vector. The caller must ensure r_hat is a unit vector.